From 0c99925f5c5f3ded5f83a5db3e265e6e01b20a75 Mon Sep 17 00:00:00 2001 From: Guillaume Ducrozet Date: Mon, 23 Feb 2015 14:18:11 +0100 Subject: [PATCH] HOTFIX: Bug corrected with respect to initialization of 3D wavefields + minor corrections - energy computation - screen output - Max. number of probes Signed-off-by: Guillaume Ducrozet --- sources/Benchmark/check_W_DYue.f90 | 18 ++++++++-------- sources/HOS/HOS-ocean.f90 | 12 ++++++----- sources/HOS/energy_calc.f90 | 8 +++---- sources/HOS/initial_condition.f90 | 20 +++++++++--------- sources/HOS/output.f90 | 21 ++++++++++++------- sources/HOS/variables_3D.f90 | 2 +- .../PostProcessing/output_post_process.f90 | 8 ++----- sources/utilities/type.f90 | 2 +- 8 files changed, 47 insertions(+), 44 deletions(-) diff --git a/sources/Benchmark/check_W_DYue.f90 b/sources/Benchmark/check_W_DYue.f90 index afb14b3..eaae129 100644 --- a/sources/Benchmark/check_W_DYue.f90 +++ b/sources/Benchmark/check_W_DYue.f90 @@ -69,7 +69,7 @@ PROGRAM check_W_DYue ! Adaptive Time Step Runge Kutta scheme TYPE(RK_parameters) :: RK_param ! -REAL(RP) :: dt_out, dt_rk4, dt_lin, dt, dt_correc, volume, energy(4) +REAL(RP) :: dt_out, dt_rk4, dt_lin, dt, dt_correc, volume, energy(3) REAL(RP) :: time_cur, time_next, h_rk, h_loc REAL(RP) :: error_old, error2 INTEGER :: idx_max(2) @@ -385,7 +385,7 @@ PROGRAM check_W_DYue volume = REAL(a_eta(1,1),RP) energy = calc_energy(a_eta, a_phis, da_eta) ! - PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4) + PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3) ! CALL fourier_2_space(a_eta, eta) CALL fourier_2_space(a_phis,phis) @@ -493,7 +493,7 @@ PROGRAM check_W_DYue volume = REAL(a_eta(1,1),RP) energy = calc_energy(a_eta, a_phis, da_eta) ! - PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4) + PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3) ! CALL fourier_2_space(a_eta, eta) CALL fourier_2_space(a_phis,phis) @@ -1211,7 +1211,7 @@ PROGRAM check_W_DYue volume = REAL(a_eta(1,1),RP) energy = calc_energy(a_eta, a_phis, da_eta) ! - PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4) + PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3) ! CALL fourier_2_space(a_eta, eta) CALL fourier_2_space(a_phis,phis) @@ -1318,7 +1318,7 @@ PROGRAM check_W_DYue volume = REAL(a_eta(1,1),RP) energy = calc_energy(a_eta, a_phis, da_eta) ! - PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4) + PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3) ! CALL fourier_2_space(a_eta, eta) CALL fourier_2_space(a_phis,phis) @@ -2182,7 +2182,7 @@ PROGRAM check_W_DYue volume = REAL(a_eta(1,1),RP) energy = calc_energy(a_eta, a_phis, da_eta) ! - PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4) + PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3) ! CALL fourier_2_space(a_eta, eta) CALL fourier_2_space(a_phis,phis) @@ -2291,7 +2291,7 @@ PROGRAM check_W_DYue volume = REAL(a_eta(1,1),RP) energy = calc_energy(a_eta, a_phis, da_eta) ! - PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4) + PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3) ! CALL fourier_2_space(a_eta, eta) CALL fourier_2_space(a_phis,phis) @@ -2307,13 +2307,13 @@ PROGRAM check_W_DYue WRITE(77,'(A,ES9.2,A,I4,A,I4)')'ZONE T = "',time_cur,'", I=',n1o2p1,', J=',n2 DO i2 = n2o2p1+1, n2 DO i1 = 1, n1o2p1 - WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), - ky(n2-i2+2), & + WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), ky_n2(i2), & ABS(a_eta(i1,i2)), ABS(a_phis(i1,i2)) ENDDO ENDDO DO i2 = 1, n2o2p1 DO i1 = 1, n1o2p1 - WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), ky(i2), & + WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), ky_n2(i2), & ABS(a_eta(i1,i2)), ABS(a_phis(i1,i2)) ENDDO ENDDO diff --git a/sources/HOS/HOS-ocean.f90 b/sources/HOS/HOS-ocean.f90 index 810c36a..00c71bf 100644 --- a/sources/HOS/HOS-ocean.f90 +++ b/sources/HOS/HOS-ocean.f90 @@ -43,7 +43,7 @@ PROGRAM HOS_ocean REAL(RP) :: delx, dely, pioxlen, pioylen, k2, thkh INTEGER :: Mo2 INTEGER :: i1, i2, j, iCPUtime, jj -REAL(RP) :: energy(4) +REAL(RP) :: energy(3) ! REAL(RP) :: t_i, t_f, time_cur, t_tot, time_next, t_i_indiv, t_f_indiv REAL(RP) :: dt_out, dt_rk4, dt, h_rk, h_loc, dt_lin @@ -284,7 +284,8 @@ PROGRAM HOS_ocean ! some variables useful for initialization ! FIXME : make it cleaner DO i2=1,n2o2p1 - DO i1=1,n1o2p1 + theta_abs(1,i2) = 0.0_rp + DO i1=2,n1o2p1 k_abs(i1,i2)=SQRT(kx(i1)*kx(i1)+(ky_n2(i2)*ky_n2(i2))) theta_abs(i1,i2)=ATAN2(ky_n2(i2),kx(i1)) IF (theta_abs(i1,i2) .LT. 0.0_rp) THEN @@ -293,7 +294,8 @@ PROGRAM HOS_ocean ENDDO ENDDO DO i2=2,n2o2p1 - DO i1=1,n1o2p1 + theta_abs(1,n2-i2+1) = 0.0_rp + DO i1=2,n1o2p1 k_abs(i1,n2-i2+2)=SQRT(kx(i1)*kx(i1)+(ky_n2(i2)*ky_n2(i2))) theta_abs(i1,n2-i2+2)=ATAN2(-ky_n2(i2),kx(i1)) IF (theta_abs(i1,n2-i2+2) .LT. 0.0_rp) THEN @@ -503,8 +505,8 @@ PROGRAM HOS_ocean ! ! Rough estimation of peak period idx = MAXLOC(abs(a_eta)) - Tp = TWOPI/(k(idx(1),idx(2))*TANH(k(idx(1),idx(2))*depth_star)) - PRINT*,'Hs_out=',4.0_rp*SQRT(energy(4))* L_out,', T_peak=',Tp * T_out + Tp = TWOPI/omega_n2(idx(1),idx(2)) + PRINT*,'Hs_out=',4.0_rp*SQRT(energy(3))* L_out,', T_peak=',Tp * T_out PRINT*,'***************************' CALL output_time_step(i_3d=i_3d, i_a=i_a_3d, i_vol=1, i_2D=i_2d, i_max=0, & time=time_cur, N_stop=FLOOR(T_stop_star / dt_out), & diff --git a/sources/HOS/energy_calc.f90 b/sources/HOS/energy_calc.f90 index 2eaf4be..516e9f0 100755 --- a/sources/HOS/energy_calc.f90 +++ b/sources/HOS/energy_calc.f90 @@ -48,7 +48,7 @@ FUNCTION calc_energy(a_eta, a_phis, da_eta) ! needs : non-dimensional gravity, g_star ! COMPLEX(CP), DIMENSION(m1o2p1,m2), INTENT(IN) :: a_eta, a_phis, da_eta -REAL(RP), DIMENSION(4) :: calc_energy +REAL(RP), DIMENSION(3) :: calc_energy ! calc_energy(:)=0.0_rp ! Potential energy @@ -81,9 +81,9 @@ FUNCTION calc_energy(a_eta, a_phis, da_eta) ! integral CALL space_2_fourier_big(temp_R_Nd, temp_C_Nd) ! -calc_energy(3) = REAL(temp_C_Nd(1,1),RP) -calc_energy(4) = calc_energy(1) + calc_energy(2) + calc_energy(3) -calc_energy(1:4) = 0.5_rp * calc_energy(1:4) +calc_energy(2) = REAL(temp_C_Nd(1,1),RP) +calc_energy(3) = calc_energy(1) + calc_energy(2) +calc_energy(1:3) = 0.5_rp * calc_energy(1:3) ! RETURN ! diff --git a/sources/HOS/initial_condition.f90 b/sources/HOS/initial_condition.f90 index b3980fa..8e9d22c 100644 --- a/sources/HOS/initial_condition.f90 +++ b/sources/HOS/initial_condition.f90 @@ -378,15 +378,15 @@ SUBROUTINE initiate_irreg(RK_param) !------------------------------------------- ! IMPLICIT NONE - +! REAL(RP) :: sigma, theta, E, Cj, pioxlen, pioylen REAL(RP) :: test, rnd, angle, angle1, angle2 -REAL(RP), DIMENSION(m1o2p1,n2) :: phi_JONSWAP -INTEGER :: iseed,i1,i2,i_initiate_NL -REAL(RP), DIMENSION(4) :: energy -TYPE(RK_parameters), INTENT(IN) :: RK_param +INTEGER :: iseed,i1,i2,i_initiate_NL +REAL(RP), DIMENSION(3) :: energy +TYPE(RK_parameters), INTENT(IN) :: RK_param COMPLEX(CP), DIMENSION(m1o2p1,m2) :: da_eta,a_eta_temp,a_phi_temp -REAL(RP),DIMENSION(m1o2p1,m2) :: DD_WW +REAL(RP), DIMENSION(m1o2p1,n2) :: phi_JONSWAP +REAL(RP),DIMENSION(m1o2p1,m2) :: DD_WW REAL(RP),DIMENSION(ikp) :: ones_k INTEGER :: i_dir_JSWP,jj REAL(RP) :: beta_min,beta_max,frr,fp_w,s_ww,Norm_DD @@ -426,12 +426,12 @@ SUBROUTINE initiate_irreg(RK_param) ELSE DO i1=1,n1o2p1 DO i2=1,n2 - DD_WW(i1,i2) =1.0_rp/beta*(cos(TWOPI*theta_abs(i1,i2)/(4*beta)))**2.0_rp + ! With this definition of angular description, theta must be in [-pi ; pi] + DD_WW(i1,i2) = 1.0_rp/beta*(cos(TWOPI*ATAN2(ky_n2(i2),kx(i1))/(4*beta)))**2.0_rp ENDDO ENDDO ENDIF ! ----- End dir function - iseed = 4*n1o2p1*n2_all ! Cj = 3.279_rp*E_cible @@ -567,7 +567,7 @@ SUBROUTINE initiate_irreg(RK_param) CALL RK_adapt_2var_3D_in_mo_lin(0, RK_param, 0.0_rp, 0.0_rp, a_phi_temp, a_eta_temp, da_eta) energy = calc_energy(a_eta, a_phis, da_eta) WRITE(*,*)'***************' - test = energy(4) + test = energy(3) ELSE IF(i_initiate_NL == 2)THEN ! FIXME : depend on wind input/dissip ? !CALL initiate_NL_o2 @@ -580,7 +580,7 @@ SUBROUTINE initiate_irreg(RK_param) ELSE CALL RK_adapt_2var_3D_in_mo_lin(0, RK_param, 0.0_rp, 0.0_rp, a_phi_temp, a_eta_temp, da_eta) energy = calc_energy(a_eta, a_phis, da_eta) - test = energy(4) + test = energy(3) ENDIF WRITE(*,*) 'E_current=', test, 'E cible =', E_cible E = E * E_cible /test diff --git a/sources/HOS/output.f90 b/sources/HOS/output.f90 index 8c82e0e..9a8dffe 100644 --- a/sources/HOS/output.f90 +++ b/sources/HOS/output.f90 @@ -64,7 +64,7 @@ SUBROUTINE init_output(i_3d, i_a, i_vol, i_2D, i_max, i_prob, i_sw) OPEN(3,file='Results/vol_energy.dat',status='unknown') CALL write_input(3) WRITE(3,'(A)')'TITLE=" 3D volume and energy "' - WRITE(3,'(A)') 'VARIABLES="t", "volume", "potential", "flexural", "kinetic", "total", "dE/E_o","E_spec_dens"' + WRITE(3,'(A)') 'VARIABLES="t", "volume", "potential", "kinetic", "total", "dE/E_o","E_spec_dens"' ENDIF ! IF (i_2D == 1) THEN @@ -136,7 +136,7 @@ SUBROUTINE init_output(i_3d, i_a, i_vol, i_2D, i_max, i_prob, i_sw) OPEN(99,file='Results/probes.dat',status='unknown') CALL write_input(9) WRITE(99,'(A)') 'TITLE="Probes records versus time"' ! - WRITE(99,'(61A)') 'VARIABLES = "time" ', ('"p'//TRIM(int2str(i1))//'" ',i1=1,nprobes) + WRITE(99,'(101A)') 'VARIABLES = "time" ', ('"p'//TRIM(int2str(i1))//'" ',i1=1,nprobes) ENDIF ! IF (i_sw == 1) THEN @@ -176,7 +176,7 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta, INTEGER, INTENT(IN) :: i_3d, i_a, i_vol, i_2D, i_max, N_stop, i_prob COMPLEX(CP), DIMENSION(m1o2p1,m2) :: a_eta, a_phis, da_eta REAL(RP), DIMENSION(m1,m2) :: eta, phis -REAL(RP) :: time, volume, energy(4), E_0(4), E_tot +REAL(RP) :: time, volume, energy(3), E_0(3), E_tot ! Local variables INTEGER :: ii, i1, i2, it REAL(RP) :: a_1, eta_mult, min_val, max_val @@ -286,9 +286,10 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta, ELSE WRITE(2,103)'ZONE T = "',time*T_out,'", I=',n1o2p1,', J=',n2 ENDIF + ! Negatives ky first DO i2 = n2o2p1+1, n2 DO i1 = 1, n1o2p1 - WRITE(2,202) kx(i1)/L_out, - ky_n2(n2-i2+2)/L_out, abs_eta(i1,i2), abs_phis(i1,i2), log_eta(i1,i2), log_phis(i1,i2) + WRITE(2,202) kx(i1)/L_out, ky_n2(i2)/L_out, abs_eta(i1,i2), abs_phis(i1,i2), log_eta(i1,i2), log_phis(i1,i2) ENDDO ENDDO DO i2 = 1, n2o2p1 @@ -321,10 +322,10 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta, IF (i_vol == 1) THEN volume = volume * xlen_star IF (n2 /= 1) volume = volume * ylen_star - IF (ABS(E_o(4)) > tiny) THEN - WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,4), ABS(energy(4)-E_0(4))/E_0(4),E_tot + IF (ABS(E_0(3)) > tiny) THEN + WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,3), ABS(energy(3)-E_0(3))/E_0(3),E_tot ELSE - WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,4), 0.0_rp,E_tot + WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,3), 0.0_rp,E_tot ENDIF 303 FORMAT(7 (ES12.5,X),ES12.5) ENDIF @@ -374,7 +375,11 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta, ENDDO ENDDO ! For the output of probes, use Hs_real and Tp_real... - WRITE(99,'(6(ES13.5,X))') time * T_out, (eta_probe(ii)* L_out, ii=1,nprobes) + WRITE(99,'(101(ES13.5,X))') time * T_out, (eta_probe(ii)* L_out, ii=1,nprobes) + IF (nprobes.GT.100) THEN + PRINT*, 'Change writing format in probes.dat file' + STOP + ENDIF ENDIF ! IF (i_sw == 1) THEN ! time=0 has to be saved... diff --git a/sources/HOS/variables_3D.f90 b/sources/HOS/variables_3D.f90 index 17aec16..4ea84a5 100644 --- a/sources/HOS/variables_3D.f90 +++ b/sources/HOS/variables_3D.f90 @@ -126,7 +126,7 @@ MODULE variables_3D REAL(RP) :: g_star, xlen_star, ylen_star, T_stop_star, f_out_star, depth_star ! ! Volume and energy -REAL(RP) :: volume, E_o(4) +REAL(RP) :: volume, E_o(3) ! ! Output numbers INTEGER :: i_3d, i_a_3d, i_2d, i_prob diff --git a/sources/PostProcessing/output_post_process.f90 b/sources/PostProcessing/output_post_process.f90 index 93b7223..cdeb2db 100644 --- a/sources/PostProcessing/output_post_process.f90 +++ b/sources/PostProcessing/output_post_process.f90 @@ -1,11 +1,7 @@ MODULE output_post_process ! -! This module contains the input related routines -! Subroutines : read_input -! read_datum -! read_blank_line -! build_format -! error_message +! This module contains the output related routines for initialization of output files and +! writing at each time-step ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! diff --git a/sources/utilities/type.f90 b/sources/utilities/type.f90 index 87b04eb..bbb4da4 100755 --- a/sources/utilities/type.f90 +++ b/sources/utilities/type.f90 @@ -21,7 +21,7 @@ MODULE type ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! -! Definition of symboles for real types (RP) and complex ones (CP) +! Definition of symbols for real types (RP) and complex ones (CP) ! ! Real numbers are simple or double precision INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6, 37) ! REAL32