diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index 8e71ef4..9295e4a 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -57,6 +57,7 @@ module namelist_mod integer :: fire_upwinding_reinit = 4 ! "numerical scheme (space) for reinitialization PDE: 1=WENO3, 2=WENO5, 3=hybrid WENO3-ENO1, 4=hybrid WENO5-ENO1" integer :: fire_lsm_band_ngp = 4 ! "number of grid points around lfn=0 that WENO5/3 is used (ENO1 elsewhere), ! for fire_upwinding_reinit=4,5 and fire_upwinding=8,9 options" + real :: reinit_pseudot_coef = 0.01 ! Coefficient for the pseudo time integer :: fast_dist_reinit_opt = 0 ! Fast distance reinitialization method (or eikonal solver): 0) None, 1) FSM integer :: fast_dist_reinit_freq = 600 ! Number of time steps to perform a reinit with fast distance reinit method @@ -77,6 +78,7 @@ module namelist_mod real :: fmoist_dt = 600.0 ! moisture model time step [s] integer :: ideal_opt = 0 ! 0) real world, 1) ideal + integer :: devel_opt = 0 ! 0) Standard nml options, 1) reads options in the devel nml block ! Objects integer :: fuel_opt = FUEL_ANDERSON ! 1) Anderson 13 @@ -154,10 +156,15 @@ module namelist_mod ! Atm block integer :: interval_atm = -1 ! Time step [s] of the atm (or frequency to read atm data if offline) integer :: kde = 1 ! Number of atm vertical levels + + ! Devel block + integer :: check_isolated_neg_lfn = 0 ! 0) Nothing, 1) Check for isolated negative values of the level set function + integer :: output_level = 0 ! 0) Standard output, >0) Specialized output contains procedure, public :: Broadcast_nml => Broadcast_nml procedure, public :: Check_nml => Check_nml procedure, public :: Initialization => Init_namelist + procedure, public :: Init_devel_block => Init_devel_block procedure, public :: Init_fire_block => Init_fire_block procedure, public :: Init_ideal_block => Init_ideal_block procedure, public :: Init_time_block => Init_time_block @@ -184,6 +191,7 @@ subroutine Broadcast_nml (this) call Broadcast_integer (this%fire_lsm_reinit_iter) call Broadcast_integer (this%fire_upwinding_reinit) call Broadcast_integer (this%fire_lsm_band_ngp) + call Broadcast_real (this%reinit_pseudot_coef) call Broadcast_integer (this%fast_dist_reinit_opt) call Broadcast_integer (this%fast_dist_reinit_freq) call Broadcast_logical (this%fire_lsm_zcoupling) @@ -202,6 +210,7 @@ subroutine Broadcast_nml (this) call Broadcast_real (this%fuelmc_c) call Broadcast_integer (this%ideal_opt) + call Broadcast_integer (this%devel_opt) call Broadcast_integer (this%fuel_opt) call Broadcast_integer (this%ros_opt) call Broadcast_integer (this%emis_opt) @@ -298,6 +307,9 @@ subroutine Broadcast_nml (this) call Broadcast_integer (this%interval_atm) call Broadcast_integer (this%kde) + ! Devel block + call Broadcast_integer (this%check_isolated_neg_lfn) + call Broadcast_integer (this%output_level) contains subroutine Broadcast_integer (val) @@ -395,6 +407,38 @@ subroutine Init_atm_block (this, file_name) end subroutine Init_atm_block + subroutine Init_devel_block (this, file_name) + + implicit none + + class (namelist_t), intent (in out) :: this + character (len = *), intent (in) :: file_name + + integer :: check_isolated_neg_lfn, output_level + integer :: unit_nml, io_stat + character (len = :), allocatable :: msg + + namelist /devel/ check_isolated_neg_lfn, output_level + + + check_isolated_neg_lfn = this%check_isolated_neg_lfn + output_level = this%output_level + + open (newunit = unit_nml, file = trim (file_name), action = 'read', iostat = io_stat) + if (io_stat /= 0) then + msg = 'Problems opening namelist file ' // trim (file_name) + call Stop_simulation (msg) + end if + + read (unit_nml, nml = devel, iostat = io_stat) + if (io_stat /= 0) call Stop_simulation ('Problems reading namelist devel block') + close (unit_nml) + + this%check_isolated_neg_lfn = check_isolated_neg_lfn + this%output_level = output_level + + end subroutine Init_devel_block + subroutine Init_fire_block (this, file_name) implicit none @@ -403,10 +447,10 @@ subroutine Init_fire_block (this, file_name) character (len = *), intent (in) :: file_name integer :: fire_print_msg, fire_upwinding, fire_lsm_reinit_iter, fire_upwinding_reinit, fire_lsm_band_ngp, & - fast_dist_reinit_opt, fast_dist_reinit_freq, fire_viscosity_ngp, wind_vinterp_opt, hinterp_opt, ideal_opt, & + fast_dist_reinit_opt, fast_dist_reinit_freq, fire_viscosity_ngp, wind_vinterp_opt, hinterp_opt, ideal_opt, devel_opt, & fuel_opt, ros_opt, fmc_opt, emis_opt, fmoist_freq real :: fire_atm_feedback, fire_viscosity, fire_lsm_zcoupling_ref, fire_viscosity_bg, fire_viscosity_band, & - fmoist_dt, fire_wind_height, frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c + fmoist_dt, fire_wind_height, frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, reinit_pseudot_coef logical :: fire_lsm_reinit, fire_lsm_zcoupling, fmoist_run, fire_is_real_perim ! ignitions @@ -426,8 +470,8 @@ subroutine Init_fire_block (this, file_name) fast_dist_reinit_opt, fast_dist_reinit_freq, fire_lsm_reinit_iter, fire_upwinding_reinit, & fire_lsm_band_ngp, fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_viscosity_bg, fire_viscosity_band, & fire_viscosity_ngp, fmoist_run, fmoist_freq, fmoist_dt, fire_wind_height, fire_is_real_perim, & - frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, ideal_opt, fuel_opt, ros_opt, fmc_opt, emis_opt, & - wind_vinterp_opt, hinterp_opt, & + frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, ideal_opt, devel_opt, fuel_opt, ros_opt, fmc_opt, emis_opt, & + wind_vinterp_opt, hinterp_opt, reinit_pseudot_coef, & ! Ignitions fire_num_ignitions, & ! Ignition 1 @@ -459,6 +503,7 @@ subroutine Init_fire_block (this, file_name) fire_lsm_reinit_iter = this%fire_lsm_reinit_iter fire_upwinding_reinit = this%fire_upwinding_reinit fire_lsm_band_ngp = this%fire_lsm_band_ngp + reinit_pseudot_coef = this%reinit_pseudot_coef fast_dist_reinit_opt = this%fast_dist_reinit_opt fast_dist_reinit_freq = this%fast_dist_reinit_freq fire_lsm_zcoupling = this%fire_lsm_zcoupling @@ -477,6 +522,7 @@ subroutine Init_fire_block (this, file_name) fuelmc_c = this%fuelmc_c ideal_opt = this%ideal_opt + devel_opt = this%devel_opt fuel_opt = this%fuel_opt ros_opt = this%ros_opt @@ -553,6 +599,7 @@ subroutine Init_fire_block (this, file_name) this%fire_lsm_reinit_iter = fire_lsm_reinit_iter this%fire_upwinding_reinit = fire_upwinding_reinit this%fire_lsm_band_ngp = fire_lsm_band_ngp + this%reinit_pseudot_coef = reinit_pseudot_coef this%fast_dist_reinit_opt = fast_dist_reinit_opt this%fast_dist_reinit_freq = fast_dist_reinit_freq this%fire_lsm_zcoupling = fire_lsm_zcoupling @@ -571,6 +618,7 @@ subroutine Init_fire_block (this, file_name) this%fuelmc_c = fuelmc_c this%ideal_opt = ideal_opt + this%devel_opt = devel_opt this%fuel_opt = fuel_opt this%ros_opt = ros_opt @@ -786,6 +834,7 @@ subroutine Init_namelist (this, file_name) call this%Init_fire_block (file_name = trim (file_name)) call this%Init_atm_block (file_name = trim (file_name)) if (this%ideal_opt > 0) call this%Init_ideal_block (file_name = trim (file_name)) + if (this%devel_opt > 0) call this%Init_devel_block (file_name = trim (file_name)) call this%Check_nml () diff --git a/physics/fire_model_mod.F90 b/physics/fire_model_mod.F90 index 4efa9d7..0e5679f 100644 --- a/physics/fire_model_mod.F90 +++ b/physics/fire_model_mod.F90 @@ -29,7 +29,6 @@ subroutine Advance_fire_model (config_flags, grid) integer :: ij, ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme real :: tbound, time_start - integer, parameter :: CHECK_ISOLATED_NEG_LFN = 0 logical, parameter :: DEBUG_LOCAL = .false. @@ -53,7 +52,8 @@ subroutine Advance_fire_model (config_flags, grid) config_flags%fire_upwinding, config_flags%fire_viscosity, config_flags%fire_viscosity_bg, config_flags%fire_viscosity_band, & config_flags%fire_viscosity_ngp, config_flags%fire_lsm_band_ngp, tbound, grid%lfn, grid%lfn_0, grid%lfn_1, grid%lfn_2, & grid%lfn_out, grid%tign_g, grid%ros, grid%uf, grid%vf, grid%dzdxf, grid%dzdyf, grid%ros_param, grid%cart_comm, & - grid%ifps, grid%ifpe, grid%jfps, grid%jfpe) + grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, grid%grad_norm_ls, grid%grad_norm_residual_sq_sum, & + grid%grad_norm_residual_sq_sum_band, grid%grad_norm_residual_rms_band) if (DEBUG_LOCAL) call Print_message ('calling Stop_if_close_to_bdy...') !$OMP PARALLEL DO & @@ -113,7 +113,7 @@ subroutine Advance_fire_model (config_flags, grid) ifds, ifde, jfds, jfde, time_start, grid%dt, grid%dx, grid%dy, config_flags%fire_upwinding_reinit, & config_flags%fire_lsm_reinit_iter, config_flags%fire_lsm_band_ngp, grid%lfn, grid%lfn_2, grid%lfn_s0, & grid%lfn_s1, grid%lfn_s2, grid%lfn_s3, grid%lfn_out, grid%tign_g, grid%cart_comm, & - grid%ifps, grid%ifpe, grid%jfps, grid%jfpe) + grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, config_flags%reinit_pseudot_coef, grid%grad_norm_reinit) if (DEBUG_LOCAL) call Print_message ('calling Copy_lfnout_to_lfn...') !$OMP PARALLEL DO & @@ -132,7 +132,7 @@ subroutine Advance_fire_model (config_flags, grid) call Do_halo_exchange_with_corners (grid%lfn, ifms, ifme, jfms, jfme, grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, N_POINTS_IN_HALO, grid%cart_comm) #endif - if (CHECK_ISOLATED_NEG_LFN == 1) call Check_isolated_negative_lfn (grid) + if (config_flags%check_isolated_neg_lfn == 1) call Check_isolated_negative_lfn (grid) if (DEBUG_LOCAL) call Print_message ('calling Ignite_prescribed_fires...') !$OMP PARALLEL DO & diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 0f3d35e..5aa7bdf 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -18,9 +18,10 @@ module level_set_mod use state_mod, only: state_fire_t, N_POINTS_IN_HALO use ignition_line_mod, only : ignition_line_t use ros_mod, only : ros_t + use constants_mod, only : PI #ifdef DM_PARALLEL - use mpi_mod, only : Do_halo_exchange, Sum_across_mpi_tasks, Max_across_mpi_tasks + use mpi_mod, only : Do_halo_exchange, Sum_across_mpi_tasks, Max_across_mpi_tasks, Min_across_mpi_tasks #endif implicit none @@ -426,7 +427,8 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & num_tiles, i_start, i_end, j_start, j_end, ts, dt, dx, dy, fire_upwinding, fire_viscosity, & fire_viscosity_bg, fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, & tbound, lfn_in, lfn_0, lfn_1, lfn_2, lfn_out, tign, ros, uf, vf, dzdxf, dzdyf, ros_model, cart_comm, & - ifps, ifpe, jfps, jfpe) + ifps, ifpe, jfps, jfpe, grad_norm_ls, grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, & + grad_norm_residual_rms_band) ! Purpose: Advance the level set function from time ts to time ts + dt @@ -438,17 +440,18 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & real, intent(in) :: fire_viscosity, fire_viscosity_bg, fire_viscosity_band real, dimension(ifms:ifme, jfms:jfme), intent (in) :: uf, vf, dzdxf, dzdyf real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign, lfn_1, lfn_2, lfn_0 - real, dimension(ifms:ifme, jfms:jfme), intent (out) :: lfn_out, ros + real, dimension(ifms:ifme, jfms:jfme), intent (out) :: lfn_out, ros, grad_norm_ls real, intent (in) :: dx, dy, ts, dt - real, intent (out) :: tbound + real, intent (out) :: tbound, grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band class (ros_t), intent (in) :: ros_model ! to store tendency (rhs of the level set pde) real, dimension(ifms:ifme, jfms:jfme) :: tend - real :: tbound2, tbound3, tbound_thread, tbound_min - integer :: i, j, ij, ifts, ifte, jfts, jfte + real :: tbound2, tbound3, tbound_thread, tbound_min, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, & + np_band, gmin, gsum + integer :: i, j, ij, ifts, ifte, jfts, jfte, np_band_local character (len = 128) :: msg - logical, parameter :: DEBUG_LOCAL = .false. + logical, parameter :: DEBUG_LOCAL = .false., PRINT_ERRORS = .false. if (DEBUG_LOCAL) call Print_message ('Entering sub Prop_level_set...') @@ -475,9 +478,15 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & ! Runge-Kutta step 1 if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 1...') - tbound_min = huge(tbound_min) + tbound_min = huge (tbound_min) + grad_norm_residual_sq_sum = 0.0 + grad_norm_residual_sq_sum_band = 0.0 + np_band = 0.0 !$OMP PARALLEL DO & - !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread, grad_norm_residual_sq_sum_local, & + !$OMP grad_norm_residual_sq_sum_band_local, np_band_local) & + !$OMP REDUCTION(MIN: tbound_min) & + !$OMP REDUCTION(+: grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, np_band) do ij = 1, num_tiles ifts = i_start(ij) ifte = i_end(ij) @@ -487,12 +496,44 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & ifms, ifme, jfms, jfme, ts, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) tbound_min = min(tbound_min, tbound_thread) + grad_norm_residual_sq_sum = grad_norm_residual_sq_sum + grad_norm_residual_sq_sum_local + grad_norm_residual_sq_sum_band = grad_norm_residual_sq_sum_band + grad_norm_residual_sq_sum_band_local + np_band = np_band + real (np_band_local) end do !$OMP END PARALLEL DO + +#ifdef DM_PARALLEL + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum, cart_comm, gsum) + grad_norm_residual_sq_sum = gsum + + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum_band, cart_comm, gsum) + grad_norm_residual_sq_sum_band = gsum + + call Sum_across_mpi_tasks (np_band, cart_comm, gsum) + np_band = gsum + + call Min_across_mpi_tasks (tbound_min, cart_comm, gmin) + tbound_min = gmin +#endif + tbound = tbound_min + grad_norm_residual_sq_sum = sqrt (grad_norm_residual_sq_sum * dx * dy) + grad_norm_residual_sq_sum_band = sqrt (grad_norm_residual_sq_sum_band * dx * dy) + if (np_band > 0.0) then + grad_norm_residual_rms_band = sqrt (grad_norm_residual_sq_sum_band / np_band) + else + grad_norm_residual_rms_band = 0.0 + end if + + if (PRINT_ERRORS) then + write (msg, '(a, 3f12.6)') 'grad_norm_residual_sq_sum rk1, band, rms_band = ', grad_norm_residual_sq_sum, & + grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band + call Print_message (msg) + end if !$OMP PARALLEL DO & !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) @@ -518,9 +559,15 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & ! Runge-Kutta step 2 if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 2...') - tbound_min = huge(tbound_min) + tbound_min = huge (tbound_min) + grad_norm_residual_sq_sum = 0.0 + grad_norm_residual_sq_sum_band = 0.0 + np_band = 0.0 !$OMP PARALLEL DO & - !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread, grad_norm_residual_sq_sum_local, & + !$OMP grad_norm_residual_sq_sum_band_local, np_band_local) & + !$OMP REDUCTION(MIN: tbound_min) & + !$OMP REDUCTION(+: grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, np_band) do ij = 1, num_tiles ifts = i_start(ij) ifte = i_end(ij) @@ -530,12 +577,44 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & ifms,ifme,jfms,jfme, ts + dt, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) - tbound_min = min(tbound_min, tbound_thread) + tbound_min = min(tbound_min, tbound_thread) + grad_norm_residual_sq_sum = grad_norm_residual_sq_sum + grad_norm_residual_sq_sum_local + grad_norm_residual_sq_sum_band = grad_norm_residual_sq_sum_band + grad_norm_residual_sq_sum_band_local + np_band = np_band + real (np_band_local) end do !$OMP END PARALLEL DO + +#ifdef DM_PARALLEL + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum, cart_comm, gsum) + grad_norm_residual_sq_sum = gsum + + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum_band, cart_comm, gsum) + grad_norm_residual_sq_sum_band = gsum + + call Sum_across_mpi_tasks (np_band, cart_comm, gsum) + np_band = gsum + + call Min_across_mpi_tasks (tbound_min, cart_comm, gmin) + tbound_min = gmin +#endif + tbound2 = tbound_min + grad_norm_residual_sq_sum = sqrt (grad_norm_residual_sq_sum * dx * dy) + grad_norm_residual_sq_sum_band = sqrt (grad_norm_residual_sq_sum_band * dx * dy) + if (np_band > 0.0) then + grad_norm_residual_rms_band = sqrt (grad_norm_residual_sq_sum_band / np_band) + else + grad_norm_residual_rms_band = 0.0 + end if + + if (PRINT_ERRORS) then + write (msg, '(a, 3f12.6)') 'grad_norm_residual_sq_sum rk2, band, rms_band = ', grad_norm_residual_sq_sum, & + grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band + call Print_message (msg) + end if !$OMP PARALLEL DO & !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) @@ -560,9 +639,15 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & ! Runge-Kutta step 3 if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 3...') - tbound_min = huge(tbound_min) + tbound_min = huge (tbound_min) + grad_norm_residual_sq_sum = 0.0 + grad_norm_residual_sq_sum_band = 0.0 + np_band = 0.0 !$OMP PARALLEL DO & - !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread, grad_norm_residual_sq_sum_local, & + !$OMP grad_norm_residual_sq_sum_band_local, np_band_local) & + !$OMP REDUCTION(MIN: tbound_min) & + !$OMP REDUCTION(+: grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, np_band) do ij = 1, num_tiles ifts = i_start(ij) ifte = i_end(ij) @@ -572,12 +657,43 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds,ifde,jfds,jfde, ifts, ifte, jfts, jfte, & ifms, ifme, jfms, jfme, ts + dt, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) tbound_min = min(tbound_min, tbound_thread) + grad_norm_residual_sq_sum = grad_norm_residual_sq_sum + grad_norm_residual_sq_sum_local + grad_norm_residual_sq_sum_band = grad_norm_residual_sq_sum_band + grad_norm_residual_sq_sum_band_local + np_band = np_band + real (np_band_local) end do !$OMP END PARALLEL DO + +#ifdef DM_PARALLEL + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum, cart_comm, gsum) + grad_norm_residual_sq_sum = gsum + + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum_band, cart_comm, gsum) + grad_norm_residual_sq_sum_band = gsum + + call Sum_across_mpi_tasks (np_band, cart_comm, gsum) + np_band = gsum + + call Min_across_mpi_tasks (tbound_min, cart_comm, gmin) + tbound_min = gmin +#endif + tbound3 = tbound_min + grad_norm_residual_sq_sum = sqrt (grad_norm_residual_sq_sum * dx * dy) + grad_norm_residual_sq_sum_band = sqrt (grad_norm_residual_sq_sum_band * dx * dy) + if (np_band > 0.0) then + grad_norm_residual_rms_band = sqrt (grad_norm_residual_sq_sum_band / np_band) + else + grad_norm_residual_rms_band = 0.0 + end if + if (PRINT_ERRORS) then + write (msg, '(a, 3f12.6)') 'grad_norm_residual_sq_sum rk3, band, rms_band = ', grad_norm_residual_sq_sum, & + grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band + call Print_message (msg) + end if !$OMP PARALLEL DO & !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) @@ -617,7 +733,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if ifds, ifde, jfds, jfde, ts, dt, dx, dy, fire_upwinding_reinit, & fire_lsm_reinit_iter, fire_lsm_band_ngp, lfn_in, lfn_2, lfn_s0, & lfn_s1, lfn_s2, lfn_s3, lfn_out, tign, cart_comm, & - ifps, ifpe, jfps, jfpe) + ifps, ifpe, jfps, jfpe, reinit_pseudot_coef, grad_norm_reinit) ! Purpose: Level-set function reinitialization ! @@ -639,7 +755,8 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_2, lfn_s0, lfn_s1, lfn_s2, lfn_s3 real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_out - real, intent (in) :: dx, dy, ts, dt + real, dimension (ifms:ifme, jfms:jfme), intent (out) :: grad_norm_reinit + real, intent (in) :: reinit_pseudot_coef, dx, dy, ts, dt real :: dt_s, threshold_hlu integer :: nts, i, j, ij, ifts, ifte, jfts, jfte @@ -682,8 +799,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if end do !$OMP END PARALLEL DO - dt_s = 0.01 * dx - dt_s = 0.0001 * dx + dt_s = reinit_pseudot_coef * dx ! iterate to solve to steady state reinit PDE ! 1 iter each time step is enoguh @@ -700,7 +816,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & lfn_s0, lfn_s3, lfn_s3, lfn_s1, 1.0 / 3.0, & ! sign funcition, initial ls, current stage ls, next stage advanced ls, RK coefficient - fire_upwinding_reinit) + fire_upwinding_reinit, grad_norm_reinit) end do !$OMP END PARALLEL DO @@ -733,7 +849,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & lfn_s0, lfn_s3, lfn_s1, lfn_s2, 1.0 / 2.0, & - fire_upwinding_reinit) + fire_upwinding_reinit, grad_norm_reinit) end do !$OMP END PARALLEL DO @@ -766,7 +882,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & lfn_s0, lfn_s3, lfn_s2, lfn_s3, 1.0, & - fire_upwinding_reinit) + fire_upwinding_reinit, grad_norm_reinit) end do !$OMP END PARALLEL DO @@ -811,7 +927,7 @@ end subroutine Reinit_level_set subroutine Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, lfn_s0, & - lfn_ini, lfn_curr, lfn_fin, rk_coeff, fire_upwinding_reinit) + lfn_ini, lfn_curr, lfn_fin, rk_coeff, fire_upwinding_reinit, grad_norm_reinit) ! Calculates right-hand-side forcing and advances a RK-stage the level-set reinitialization PDE @@ -822,6 +938,7 @@ subroutine Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & integer, intent (in) :: fire_upwinding_reinit real, dimension (ifms:ifme, jfms:jfme), intent (in) :: lfn_s0, lfn_ini, lfn_curr real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_fin + real, dimension (ifms:ifme, jfms:jfme), intent (out) :: grad_norm_reinit real, intent (in) :: dx, dy, dt_s, threshold_hlu, rk_coeff integer :: i, j @@ -930,6 +1047,7 @@ subroutine Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & end select end if grad = sqrt (diff2x * diff2x + diff2y * diff2y) + grad_norm_reinit(i, j) = grad tend_r = lfn_s0(i, j) * (1.0 - grad) lfn_fin(i, j) = lfn_ini(i, j) + (dt_s * rk_coeff) * tend_r end do @@ -969,7 +1087,8 @@ end subroutine Update_ignition_times subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfms, jfme, & t, dt, dx, dy, fire_upwinding, fire_viscosity, fire_viscosity_bg, & - fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, lfn, tbound, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, lfn, tbound, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) ! compute the right hand side of the level set equation @@ -977,16 +1096,18 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm integer, intent (in) :: ifms, ifme, jfms, jfme, its, ite, jts, jte, ids, ide, jds, jde, & fire_upwinding,fire_viscosity_ngp, fire_lsm_band_ngp + integer, intent (out) :: np_band_local real, intent (in) :: fire_viscosity, fire_viscosity_bg, fire_viscosity_band, t, dt, dx, dy real, dimension(ifms:ifme, jfms:jfme), intent (in) :: uf, vf, dzdxf, dzdyf real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: lfn - real, dimension(ifms:ifme, jfms:jfme), intent (out) :: tend, ros - real, intent (out) :: tbound + real, dimension(ifms:ifme, jfms:jfme), intent (out) :: tend, ros, grad_norm_ls + real, intent (out) :: tbound, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local class (ros_t), intent (in) :: ros_model real, parameter :: EPS = epsilon (0.0), TOL = 100.0 * EPS real :: difflx, diffly, diffrx, diffry, diffcx, diffcy, & - diff2x, diff2y, grad, & + diff2x, diff2y, grad, mask, diff2x_eno, diff2y_eno, & + diff2x_weno, diff2y_weno, transition_width, band_width, & scale, nvx, nvy, a_valor, signo_x, signo_y, threshold_hll, & threshold_hlu, threshold_av, fire_viscosity_var integer :: i, j @@ -1006,6 +1127,9 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm if (DEBUG_LOCAL) call Print_message ('starting ij loops') tbound = 0.0 + grad_norm_residual_sq_sum_local = 0.0 + grad_norm_residual_sq_sum_band_local = 0.0 + np_band_local = 0 do j = jts, jte do i = its, ite ! one sided differences @@ -1013,51 +1137,59 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm difflx = (lfn(i, j) - lfn(i - 1, j)) / dx diffry = (lfn(i, j + 1) - lfn(i, j)) / dy diffly = (lfn(i, j) - lfn(i, j - 1)) / dy - ! twice central difference - diffcx = difflx + diffrx - diffcy = diffly + diffry + + ! central differences + diffcx = 0.5 * (difflx + diffrx) + diffcy = 0.5 * (diffly + diffry) + ! use eno1 near domain boundaries if (i < ids + BDY_ENO1 .or. i > ide - BDY_ENO1 .or. & j < jds + BDY_ENO1 .or. j > jde - BDY_ENO1) then diff2x = Select_eno (difflx, diffrx) diff2y = Select_eno (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) + else select case (fire_upwinding) - ! none case (0) - grad = sqrt (diffcx ** 2 + diffcy ** 2) + ! none + diff2x = diffcx + diff2y = diffcy + grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! standard case (1) + ! First-order sign-consistent upwind scheme (component wise). Standard diff2x = Select_upwind (difflx, diffrx) diff2y = Select_upwind (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! godunov per osher/fedkiw case (2) + ! Component-wise Godunov upwind (1st order) derivative selection diff2x = Select_godunov (difflx, diffrx) diff2y = Select_godunov (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! ENO1 case (3) + ! First order ENO (minmod limited) componenet-wise gradient diff2x = Select_eno (difflx, diffrx) diff2y = Select_eno (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! Sethian - twice stronger pushdown of bumps - case(4) - grad = sqrt (max (difflx, 0.0) ** 2 + min (diffrx, 0.0) ** 2 & - + max (diffly, 0.0) ** 2 + min(diffry, 0.0) ** 2) - ! 2nd order - case(5) + case (4) + ! Classical Godunov discretization of the Hamiltonian + grad = sqrt (max (difflx, 0.0) ** 2 + min (diffrx, 0.0) ** 2 + & + max (diffly, 0.0) ** 2 + min (diffry, 0.0) ** 2) + diff2x = max (difflx, 0.0) - min (diffrx, 0.0) + diff2y = max (diffly, 0.0) - min (diffry, 0.0) + + case (5) + ! 2nd order central differences diff2x = Select_2nd (dx, lfn(i, j), lfn(i - 1, j), lfn(i + 1, j)) diff2y = Select_2nd (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j + 1)) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! WENO3 - case(6) + case (6) + ! 3rd order WENO a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), & @@ -1070,8 +1202,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm lfn(i, j + 1), lfn(i, j + 2), signo_y) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! WENO5 - case(7) + case (7) + ! 5th order WENO a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j)+ & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) @@ -1082,8 +1214,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm lfn(i ,j - 3), lfn(i, j + 1), lfn(i, j + 2), lfn(i, j + 3), signo_y) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! WENO3/ENO1 - case(8) + case (8) + ! WENO3/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) @@ -1102,8 +1234,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) end if - ! WENO5/ENO1 - case(9) + case (9) + ! WENO5/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & Select_4th (dy,lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) @@ -1122,6 +1254,43 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) end if + case (10) + ! WENO5/ENO1 + blending zone + band_width = 2 * threshold_hlu + transition_width = threshold_hlu + if (abs(lfn(i, j)) <= band_width) then + if (abs(lfn(i, j)) <= band_width - transition_width) then + mask = 1.0 + else + mask = 0.5 * (1.0 + cos (PI * (abs (lfn(i, j)) -(band_width - transition_width)) / transition_width)) + end if + else + mask = 0.0 + end if + + diff2x_eno = Select_eno (difflx, diffrx) + diff2y_eno = Select_eno (diffly, diffry) + + if (mask > 0.0) then + a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & + Select_4th (dy,lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) + signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), & + lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) + signo_y = a_valor * Select_4th (dy, lfn(i, j), lfn(i, j - 1), & + lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) + diff2x_weno = Select_weno5 (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), & + lfn(i - 3, j), lfn(i + 1, j), lfn(i + 2, j), lfn(i + 3, j), signo_x) + diff2y_weno = Select_weno5 (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), & + lfn(i, j - 3), lfn(i, j + 1), lfn(i, j + 2), lfn(i, j + 3), signo_y) + diff2x = mask * diff2x_weno + (1.0 - mask) * diff2x_eno + diff2y = mask * diff2y_weno + (1.0 - mask) * diff2y_eno + else + diff2x = diff2x_eno + diff2y = diff2y_eno + end if + + grad = sqrt (diff2x * diff2x + diff2y * diff2y) + case default !$omp critical write (msg, '(a, i2)') 'Unknown upwinding option in level set : ', fire_upwinding @@ -1131,6 +1300,13 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm end select end if + grad_norm_ls(i, j) = grad + grad_norm_residual_sq_sum_local = grad_norm_residual_sq_sum_local + (grad - 1.0) ** 2 + if (abs(lfn(i, j)) < threshold_hlu) then + grad_norm_residual_sq_sum_band_local = grad_norm_residual_sq_sum_band_local + (grad - 1.0) ** 2 + np_band_local = np_band_local + 1 + end if + ! Calc normal scale = sqrt (grad ** 2.0 + EPS) nvx = diff2x / scale diff --git a/share/mpi_mod.F90 b/share/mpi_mod.F90 index 27aa26e..b8a1280 100644 --- a/share/mpi_mod.F90 +++ b/share/mpi_mod.F90 @@ -7,7 +7,7 @@ module mpi_mod private public :: Calc_tasks_in_x_and_y, Calc_patch_dims, Gather_var2d, Do_halo_exchange, Do_halo_exchange_with_corners, & - Max_across_mpi_tasks, Sum_across_mpi_tasks, Distribute_var2d + Max_across_mpi_tasks, Sum_across_mpi_tasks, Distribute_var2d, Min_across_mpi_tasks contains @@ -538,6 +538,29 @@ subroutine Max_across_mpi_tasks (local_max, cart_comm, global_max) end subroutine Max_across_mpi_tasks + subroutine Min_across_mpi_tasks (local_min, cart_comm, global_min) + +#ifdef DM_PARALLEL + use mpi +#endif + + implicit none + + real, intent (in) :: local_min + integer, intent (in) :: cart_comm + real, intent (out) :: global_min + + integer :: ierr + + +#ifdef DM_PARALLEL + call MPI_Allreduce (local_min, global_min, 1, MPI_REAL, MPI_MIN, cart_comm, ierr) +#else + global_min = local_min +#endif + + end subroutine Min_across_mpi_tasks + subroutine Sum_across_mpi_tasks (local_sum, cart_comm, global_sum) #ifdef DM_PARALLEL diff --git a/state/state_mod.F90 b/state/state_mod.F90 index c0f5e61..6f574f8 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -73,6 +73,8 @@ module state_mod real, dimension(:, :), allocatable :: nfuel_cat ! "fuel data" real, dimension(:, :), allocatable :: fuel_time ! "fuel" real, dimension(:, :), allocatable :: emis_smoke + real, dimension(:, :), allocatable :: grad_norm_ls ! Gracient norm of the level set function used to propagate level set function + real, dimension(:, :), allocatable :: grad_norm_reinit ! Gracient norm of the level set function used to reinitilize the level set function class (fuel_t), allocatable :: fuels class (ros_t), allocatable :: ros_param @@ -98,6 +100,14 @@ module state_mod integer :: ny ! "number of latitudinal grid points" "1" real :: cen_lat, cen_lon + ! Performance stats + real :: grad_norm_residual_sq_sum + real :: grad_norm_residual_sq_sum_band + real :: grad_norm_residual_rms_band + + ! Output + integer :: output_level + ! For MPI tasks integer :: cart_comm ! The MPI communicator with the domain decomposition integer :: ntasks ! Number of MPI tasks @@ -173,6 +183,8 @@ subroutine Allocate_vars (this, ifms, ifme, jfms, jfme) allocate (this%dzdyf(ifms:ifme, jfms:jfme)) allocate (this%nfuel_cat(ifms:ifme, jfms:jfme)) allocate (this%emis_smoke(ifms:ifme, jfms:jfme)) + allocate (this%grad_norm_ls(ifms:ifme, jfms:jfme)) + allocate (this%grad_norm_reinit(ifms:ifme, jfms:jfme)) end subroutine Allocate_vars @@ -594,6 +606,9 @@ subroutine Init_domain (this, config_flags, geogrid, & if (DEBUG_LOCAL) call Print_message (' Setting clock...') call this%Set_time_stamps (config_flags) + ! Output + this%output_level = config_flags%output_level + if (DEBUG_LOCAL) call this%Print() if (DEBUG_LOCAL) call Print_message ('Leaving Init_domain...') @@ -944,6 +959,14 @@ subroutine Save_state (this) call Add_netcdf_var_mpi (file_output, this%nx, this%ny, this%ifps, this%ifpe, this%jfps, this%jfpe, 'nfuel_cat', & this%nfuel_cat(this%ifps:this%ifpe, this%jfps:this%jfpe)) + if (this%output_level > 0) then + call Add_netcdf_var_mpi (file_output, this%nx, this%ny, this%ifps, this%ifpe, this%jfps, this%jfpe, 'grad_norm_ls', & + this%grad_norm_ls(this%ifps:this%ifpe, this%jfps:this%jfpe)) + + call Add_netcdf_var_mpi (file_output, this%nx, this%ny, this%ifps, this%ifpe, this%jfps, this%jfpe, 'grad_norm_reinit', & + this%grad_norm_reinit(this%ifps:this%ifpe, this%jfps:this%jfpe)) + end if + if (DEBUG_LOCAL) call Print_message ('Leaving Save_state...') end subroutine Save_state diff --git a/tests/test7/namelist.fire b/tests/test7/namelist.fire index d6931c7..7c02c76 100644 --- a/tests/test7/namelist.fire +++ b/tests/test7/namelist.fire @@ -37,5 +37,6 @@ fire_viscosity = 0.4, ! artificial viscosity in level set method (max 1, needed with fire_upwinding=0) fire_upwinding = 9, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian wind_vinterp_opt = 0, + reinit_pseudot_coef = 0.0001 / diff --git a/tests/test8/namelist.fire b/tests/test8/namelist.fire index 502787d..62d4b6f 100644 --- a/tests/test8/namelist.fire +++ b/tests/test8/namelist.fire @@ -39,5 +39,6 @@ fmoist_run = .true., ! Run fuel moisture model fmoist_freq = 1, ! Fuel moisture model updated every time step wind_vinterp_opt = 0 + reinit_pseudot_coef = 0.0001 / diff --git a/tests/testx/namelist.fire b/tests/testx/namelist.fire index 5dfd895..0346d66 100644 --- a/tests/testx/namelist.fire +++ b/tests/testx/namelist.fire @@ -37,6 +37,7 @@ fire_viscosity = 0.4, ! artificial viscosity in level set method (max 1, needed with fire_upwinding=0) fire_upwinding = 9, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian wind_vinterp_opt = 0 + reinit_pseudot_coef = 0.0001 / &ideal