Skip to content

Commit

Permalink
correct bugs for GQM to run correctly
Browse files Browse the repository at this point in the history
  • Loading branch information
@mickaelaccensi
mickaelaccensi committed Nov 21, 2023
1 parent 1f928aa commit 8ff8d22
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 23 deletions.
29 changes: 15 additions & 14 deletions manual/eqs/NL1.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,37 +55,38 @@ \subsubsection{~$S_{nl}$: Discrete Interaction Approximation (\dia)} \label{sec:
\sin(\delta_{\theta,3})&=&\sin(\delta_{\theta,2}) (1-\lambda)^2/(1+\lambda)^2.
\end{eqnarray}

For these quadruplets, each source term value
$S_{nl}(\bk)$ corresponding to each discrete $(f_r,\theta)$
we compute the three contributions that correspond to the situation in which $\bk$ takes the role of $\bk$,$\bk_{2,+}$, $\bk_{2,-}$, $\bk_{3,+}$ and $\bk_{3,-}$ in the quadruplet, namely the full source term is
Hence for any $\bk$ one quadruplet selects $\bk_{2,+}$ and $\bk_{3,+}$, and the other quadruplet selects its mirror image
$\bk_{2,-}$, $\bk_{2,-}$. Because there are 3 different components interacting in the two DIA-selected quadruplets, any discrete spectral component $(f_r,\theta)$ is actually involved in 6 quadruplets and directly exchanges energy with 12 other components $(f_r',\theta')$. Because the values of $f'_r$ and $\theta'$ do not fall exacly on other discrete components, the spectral density is interpolated using a bilinear interpolation, so that each source term value
$S_{nl}(\bk)$ contains the direct exchange of energy with 48 other discrete components.
we compute the three contributions that correspond to the situation in which $\bk$ takes the role of $\bk$,$\bk_{2,+}$, $\bk_{2,-}$, $\bk_{3,+}$ and $\bk_{3,-}$ in the quadruplet, namely the full source term is, without making explicit that bilinear interpolation,
\begin{eqnarray}
S_{\mathrm{nl}}(\bk) &=& -2 \left[\delta S_{\mathrm{nl}}(\bk,\bk_2,\bk_3,+)+\delta S_{\mathrm{nl}}(\bk,\bk_2,\bk_3,-)\right] \nonumber \\
& & + \delta S_{\mathrm{nl}}(\bk_4,\bk,\bk_5,+) + \delta S_{\mathrm{nl}}(\bk_6,\bk,\bk_7,-) \\
& & + \delta S_{\mathrm{nl}}(\bk_8,\bk_9,\bk, +) + \delta S_{\mathrm{nl}}(\bk_{10},\bk_{11},\bk, -) . \label{eq:diasum}
S_{\mathrm{nl}}(\bk) &=& -2 \left[\delta S_{\mathrm{nl}}(\bk,\bk_{2,+},\bk_{3,+})+\delta S_{\mathrm{nl}}(\bk,\bk_{2,-},\bk_{3,-})\right] \nonumber \\
& & + \delta S_{\mathrm{nl}}(\bk_4,\bk,\bk_5) + \delta S_{\mathrm{nl}}(\bk_6,\bk,\bk_7) \\
& & + \delta S_{\mathrm{nl}}(\bk_8,\bk_9,\bk) + \delta S_{\mathrm{nl}}(\bk_{10},\bk_{11},\bk) . \label{eq:diasum}
\end{eqnarray}
with elementary contributions given by
where the geometry of the quadruplet $(\bk_4,\bk_4,\bk,\bk_5)$ is obtained from that of $(\bk,\bk,\bk_{2,+},\bk_{3,+})$ by a dilation by a factor $(1+\lambda)^2$ and rotation by the angle $\delta_{\theta,2}$; $(\bk_6,\bk_6,\bk,\bk_7)$ has the same dilation but the opposite rotation; $(\bk_8,\bk_8,\bk_9,\bk)$ is dilated by a factor $(1-\lambda)^2$ and rotated by the angle $-\delta_{\theta,3}$: and $(\bk_{10},\bk_{10},\bk_{11},\bk)$ is dilated by the same factor and rotated by the opposite angle.


The elementary contributions $\delta S_{\mathrm{nl}}(\bk_l,\bk_m,\bk_n)$ are given by
%----------------------------%
% Nonlinear interactions DIA %
%----------------------------%
% eq:snl_dia

\begin{equation}
\delta S_{\mathrm{nl}}(\bk,\bk_2,\bk_3,s) = \frac{C}{g^4} f_{r,1}^{11} \left [ F^2 \left ( \frac{F_{2,s}}{(1+\lambda_{nl})^4} +
\frac{F_{3,s}}{(1-\lambda_{nl})^4} \right ) - \frac{2 F F_{2,s} F_{3,s}}{(1-\lambda_{nl}^2)^4} \right] ,
\delta S_{\mathrm{nl}}(\bk_l,\bk_m,\bk_n) = \frac{C}{g^4} f_{r,l}^{11} \left [ F_l^2 \left ( \frac{F_m}{(1+\lambda)^4} +
\frac{F_n}{(1-\lambda)^4} \right ) - \frac{2 F_l F_m F_n}{(1-\lambda^2)^4} \right] ,
\label{eq:snl_dia}
\end{equation}
where $s=+$ or $s=-$ is a sign index, and the spectral densities are $F = F(f_{r} ,\theta)$, $F_{2,+} = F(f_{r,2} ,\theta + \delta_{\theta,2})$, $F_{2,-} = F(f_{r,2} ,\theta - \delta_{\theta,2})$, etc.
where the spectral densities are $F_l = F(f_{r,l} ,\theta_l)$, etc.
$C$ is a proportionality constant that was tuned to reproduce the inverse energy cascade. Default values for different source term packages are presented in Table~\ref{tab:snl_par}.
As a result, when accounting for the two quadruplet configurations, the source term at $\bk$ includes the interactions with
10 other spectral components. Besides, because $f_{r,2}$ and $f_{r,3}$ nor $\theta_{2,\pm} $ and $\theta_{3,\pm} $ fall on discretized frequencies and directions, the spectral densities are bilinearly interpolated, which involves 4 discrete spectral components for each of these 10 components.



% tab:snl_par

\begin{table} \begin{center}
\begin{tabular}{|l|c|c|} \hline
& $\lambda_{nl}$ & $C$ \\ \hline
& $\lambda$ & $C$ \\ \hline
ST6 & 0.25 & $3.00 \; 10^7$ \\ \hline
\wam-3 & 0.25 & $2.78 \; 10^7$ \\ \hline
ST4 (Ardhuin et al.)& 0.25 & $2.50 \; 10^7$ \\ \hline
Expand Down
11 changes: 11 additions & 0 deletions manual/manual.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3665,6 +3665,17 @@ @article{art:DC23
year = {2023}
}

@ARTICLE{Webb1978,
author = "D. J. Webb",
title = "Nonlinear transfer between sea waves",
journal = DSR,
volume = 25,
pages = "279--298",
year = 1978,
where="paper",
}


@ARTICLE{Lavrenov2001,
author = "Igor V. Lavrenov",
title = "Effect of wind wave parameter fluctuation on the nonlinear spectrum evolution",
Expand Down
6 changes: 5 additions & 1 deletion manual/sys/files_w3.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -506,11 +506,15 @@ \subsubsection{~Wave model modules} \label{sec:wave_mod}
\end{flist}

\noindent
Nonlinear interaction module (\dia) \hfill {\file w3snl1md.ftn}
Nonlinear interaction module (\dia or GQM) \hfill {\file w3snl1md.ftn}

\begin{flisti}
\fit{w3snl1}{Calculation of $S_{nl}$.}
\fit{insnl1}{Initialization for $S_{nl}$.}
\fit{w3snlgqm}{Calculation of $S_{nl}$.}
\fit{w3scouple}{Calculation of coupling coefficient.}
\fit{gauleg}{Calculation of Gauss-Legendre quadrature coefficients.}
\fit{INSNLGQM}{Initialization for $S_{nl}$ with GQ method.}
\end{flisti}

\noindent
Expand Down
22 changes: 14 additions & 8 deletions model/src/w3snl1md.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -825,6 +825,8 @@ SUBROUTINE W3SNLGQM(A,CG,WN,DEPTH,TSTOTn,TSDERn)
USE CONSTANTS, ONLY: TPI
USE W3GDATMD, ONLY: SIG, NK , NTH , DTH, XFR, FR1, GQTHRSAT, GQAMP

IMPLICIT NONE

REAL, intent(in) :: A(NTH,NK), CG(NK), WN(NK)
REAL, intent(in) :: DEPTH
REAL, intent(out) :: TSTOTn(NTH,NK), TSDERn(NTH,NK)
Expand Down Expand Up @@ -883,8 +885,8 @@ SUBROUTINE W3SNLGQM(A,CG,WN,DEPTH,TSTOTn,TSDERn)
! Gamma_max=1.3 (JFMAX>NF) TO OBTAIN IMPROVED RESULTS
! Note by Fabrice Ardhuin: this appears to give the difference in tail benaviour with Gerbrant's WRT
!=======================================================================
JFMIN= 1-INT(LOG(1.0D0)/LOG(RAISF))
JFMAX=NF+INT(LOG(1.3D0)/LOG(RAISF))
JFMIN=MAX(1-INT(LOG(1.0D0)/LOG(RAISF)),1)
JFMAX=MIN(NF+INT(LOG(1.3D0)/LOG(RAISF)),NK)
!
!=======================================================================
! COMPUTES THE SPECTRUM THRESHOLD VALUES (BELOW WHICH QNL4 IS NOT
Expand Down Expand Up @@ -1065,7 +1067,7 @@ SUBROUTINE W3SNLGQM(A,CG,WN,DEPTH,TSTOTn,TSDERn)
TEMP=(TB_TPM(IQ_OM2,JT1,JF1)*(( F(JT1P2P,JFM2)*CF2 *F(JT1P3M,JFM3)*CF3)* &
(F(JT,JFM0 )*CF0*TB_V14(JF1)+F(JT1P ,JFM1)*CF1) &
-SP0*SP1P*(SP1P2P*V3_4+SP1P3M*V2_4))+T_2M3P*(AUX05*AUX01-AUX02*AUX06)) *CP0
WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT, F(JT,JFM0)
WRITE(995,'(3I3,3E12.3)') ICONF,JF,JT, F(JT,JFM0)
TEMP=(Q_2P3M+Q_2M3P) *CP1
WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT,JT1P, JFM1,AUX00 *CP1, F(JT1P,JFM1),TSTOT(JT1P,JFM1)
WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT,JT1P2P,JFM2,-Q_2P3M*CP2,F(JT1P2P,JFM2),TSTOT(JT1P2P,JFM2)
Expand Down Expand Up @@ -1219,6 +1221,8 @@ FUNCTION COUPLE(XK1 ,YK1 ,XK2 ,YK2 ,XK3 ,YK3 ,XK4 ,YK4)
!/ ------------------------------------------------------------------- /
USE CONSTANTS, ONLY: GRAV
!
IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: XK1 , YK1 , XK2 , YK2
DOUBLE PRECISION, INTENT(IN) :: XK3 , YK3
DOUBLE PRECISION, INTENT(IN) :: XK4 , YK4
Expand Down Expand Up @@ -1305,6 +1309,7 @@ SUBROUTINE GAULEG (W_LEG ,X_LEG ,NPOIN)
!/ ------------------------------------------------------------------- /
!.....VARIABLES IN ARGUMENT
! """"""""""""""""""""
IMPLICIT NONE
INTEGER , INTENT(IN) :: NPOIN
DOUBLE PRECISION ,INTENT(INOUT) :: W_LEG(NPOIN) , X_LEG(NPOIN)
!
Expand Down Expand Up @@ -1552,6 +1557,7 @@ SUBROUTINE INSNLGQM
#ifdef W3_S
CALL STRACE (IENT, 'INSNLGQM')
#endif
IMPLICIT NONE
!.....LOCAL VARIABLES
INTEGER JF , JT , JF1 , JT1 , NF1P1 , IAUX , NT , NF , IK
INTEGER IQ_TE1 , IQ_OM2 , LBUF , DIMBUF , IQ_OM1 , NQ_TE1 , NCONFM
Expand Down Expand Up @@ -2084,10 +2090,7 @@ SUBROUTINE INSNLGQM
AUX=0.0D0
DO JT1=1,GQNT1
DO IQ_OM2=1,GQNQ_OM2
AAA=TB_FAC(IQ_OM2,JT1,JF1)*TB_TPM(IQ_OM2,JT1,JF1)
IF (AAA.GT.AUX) AUX=AAA
CCC=TB_FAC(IQ_OM2,JT1,JF1)*TB_TMP(IQ_OM2,JT1,JF1)
IF (CCC.GT.AUX) AUX=CCC
AUX=MAX(AUX,TB_FAC(IQ_OM2,JT1,JF1)*TB_TPM(IQ_OM2,JT1,JF1),TB_FAC(IQ_OM2,JT1,JF1)*TB_TMP(IQ_OM2,JT1,JF1))
ENDDO
ENDDO
MAXCLA(JF1)=AUX
Expand All @@ -2099,6 +2102,7 @@ SUBROUTINE INSNLGQM
DO JF1=1,GQNF1
IF (MAXCLA(JF1).GT.AUX) AUX=MAXCLA(JF1)
ENDDO

TEST1=SEUIL1*AUX
!
!.....Set to zero the coupling coefficients not used
Expand Down Expand Up @@ -2128,7 +2132,9 @@ SUBROUTINE INSNLGQM
!
!..... counts the fraction of the eliminated configurations
ELIM=(1.D0-DBLE(NCONF)/DBLE(NCONFM))*100.D0
! WRITE(994,*) 'NCONF:',NCONF,ELIM
#ifdef W3_TGQM
WRITE(994,*) 'NCONF, ELIM FRACTION:',NCONF,ELIM
#endif
END SUBROUTINE INSNLGQM
!/
!/ End of module W3SNL1MD -------------------------------------------- /
Expand Down

0 comments on commit 8ff8d22

Please sign in to comment.