TELEMAC-MASCARET Forum: Wilcock and Crowe subroutine modification

Welcome, Guest

TOPIC: Wilcock and Crowe subroutine modification - segfault

Wilcock and Crowe subroutine modification - segfault 1 year 6 months ago #42541

Renault
OFFLINE
Senior Boarder
Posts: 120
Thank you received: 33

Hi all,

I am attempting to incorporate some modifications to the BEDLOAD_WILCOCK_CROWE_GAIA subroutine as described in Recking et al. 2016. These modifications are twofold:

Modifying the exponent and intercept;
Correcting shear stress.

Item 1 is simple enough, as it is a matter of changing the coefficients in the subroutine (exponent 7.5 becomes 2.0 and intercept 1.35 becomes 1.10). Item 2 has proven to be much more challenging and gives a SIGSEGV error.
The equation to be implemented, and calculated at each node, is:
τ' = 17*(SD₆₅)^1/4*U^3/2
With τ' as corrected shear stress (Pa), S as local bed slope (m/m), D₆₅ as representative grain size (mm), and U as local velocity (m/s).

I have included the subroutine file as an attachment, but I will also include the relevant portions here.
Declaration of the name of the subroutine and variables:

!                   *************************************
                    SUBROUTINE BEDLOAD_WILCOCK_CROWE_GAIA
!                   *************************************
!
     &(TOB, MU, ACLADM, DCLA, RATIO_SAND, GRAV, XMVE, XMVS, SANFRA, QSC,
     & ACP, SLOPEFF,COEFCR, DZFDX, DZFDY, UNORM)

Parameters read for DICO (I assume):

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!>@param[in]     TOB        Bed shear stress
!>@param[in]     MU         Skin friction correction factor for bed roughness
!>@param[in]     ACLADM     Mean diameter of active-layer
!>@param[in]     DCLA       Sediment grain diameter
!>@param[in]     RATIO_SAND Ratio of sand to all sands, for isand,ilayer,ipoin
!>@param[in]     GRAV       Acceleration of gravity
!>@param[in]     XMVE       Water density
!>@param[in]     XMVS       Sand density
!>@param[in]     SANFRA     Sand fraction for Wilcock & Crowe transport formula
!>@param[in,out] QSC        Bed load transport
!>@param[in,out] ACP        Coefficient to modify Shields reference
!>@param[in]     SLOPEFF    Logical, sloping bed effect or not
!>@param[in]     DZFDX      Bottom slope in the x-direction
!>@param[in]     DZFDY      Bottom slope in the y-direction
!>@param[in]     UNORM      Norm of the mean flow velocity
!>@param[in,out] COEFPN     Correction critical Shields for sloping bed effect
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Use modules:

USE BIEF
      USE BIEF_DEF
      USE INTERFACE_GAIA,
     &      EX_BEDLOAD_WILCOCK_CROWE_GAIA => BEDLOAD_WILCOCK_CROWE_GAIA
      USE DECLARATIONS_SPECIAL
      IMPLICIT NONE

Declaring input/output variables:

INTEGER,          INTENT(IN)    :: SLOPEFF
      TYPE(BIEF_OBJ),   INTENT(INOUT) :: QSC
      DOUBLE PRECISION, INTENT(IN)    :: XMVE, XMVS, GRAV, DCLA,
     &                                   RATIO_SAND(QSC%DIM1)
      DOUBLE PRECISION, INTENT(IN)    :: SANFRA(QSC%DIM1)
      TYPE(BIEF_OBJ),   INTENT(INOUT) :: ACP ! WORK ARRAY T1
      TYPE(BIEF_OBJ),   INTENT(IN)    :: ACLADM,COEFCR
      TYPE(BIEF_OBJ),   INTENT(IN)    :: UNORM
      TYPE(BIEF_OBJ),   INTENT(IN)    :: DZFDX,DZFDY
      TYPE(BIEF_OBJ),   INTENT(IN)    :: TOB, MU

Local variables:

INTEGER I
      DOUBLE PRECISION TORM, TOBPR, TORI, TORATIO, WI, COEFB, WCC, DZF
!     TOBPR = BED SHEAR STRESS PRIME, CORR FACTOR FROM WILCOCK ETAL 2009
!     DZF = SLOPE NORM DERIVED FROM BEDLOAD_EFFPNT_GAIA

Actual "meat" of the function, starting after where TORI is set, where:

the local slope S (DZF) magnitude is calculated from local X (DZFDX%R(I)) and local Y (DZFDY%R(I)) components;
the local 65th percentile diameter D₆₅ is approximated with mean diameter (ACLADM%R(I)) and multiplied by 1000 to convert from m to mm;
the local velocity magnitude (UNORM%R(I)) is obtained;
the corrected sheaer stress (TOBPR) is calculated using slope, diameter, velocity, and coefficients/exponents.

In addition, certain parameters are changed at the end.

!     CORRECTED FORMULA FROM WILCOCK ET AL 2009 EQN 2.15
!     CORRECTED FORMULA USES D65 WHICH IS NOT REALLY CALCULATED
!     ELSEWHERE IN THE MODEL. USE MEAN DIAMETER ACLADM AS APPROX OF D65
!     FLOW VELOCITY: UNORM, SLOPE: DZF
        DZF=SQRT(DZFDX%R(I)**2+DZFDY%R(I)**2)
        TOBPR = 17.D0*((DZF*ACLADM%R(I)*1.D3)**0.25D0)*
     &  (UNORM%R(I)**1.5D0)
        TORATIO = TOBPR*MU%R(I)/TORI
!     CHANGED 1.35 TO 1.10 AND 7.5 TO 2.0 AS PER RECKING ET AL 2016
        IF (TORATIO.LT.1.10D0) THEN
          WI = 2.D-3*(TORATIO**2.0D0)
        ELSE
          WI = 14.D0*((1.D0-0.894D0/SQRT(TORATIO))**4.5D0)
        ENDIF

My suspicion is that the slope or velocity are not accessed correctly, hence the segfault, but I'm not sure how to access them otherwise. The actual error message, repeated to the number of threads I have (testing in parallel) is this:

forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image              PC                Routine            Line        Source
out_user_fortran   000000000042C3C9  for__signal_handl     Unknown  Unknown
libpthread-2.30.s  0000150C8DA710F0  Unknown               Unknown  Unknown
out_user_fortran   000000000040BA04  bedload_wilcock_c     Unknown  Unknown
libgaia4telemac2d  0000150C90139ED3  bedload_formula_g     Unknown  Unknown
libgaia4telemac2d  0000150C90197752  init_transport_ga     Unknown  Unknown
libgaia4telemac2d  0000150C901A338C  gaia_init_            Unknown  Unknown
libtelemac2d.so    0000150C9036A608  telemac2d_init_       Unknown  Unknown
libtelemac2d.so    0000150C903BA2BC  telemac2d_            Unknown  Unknown
libtelemac2d.so    0000150C9025B577  MAIN__                Unknown  Unknown
out_user_fortran   000000000040B8D2  Unknown               Unknown  Unknown
libc-2.30.so       0000150C8D77DE1B  __libc_start_main     Unknown  Unknown
out_user_fortran   000000000040B7EA  Unknown               Unknown  Unknown

Any insights into what I could have done wrong, and how to correct it, are more than welcome! As may be evident, I'm still new to FORTRAN, let alone the TELEMAC way of using it.

Many thanks in advance for any insights,
André Renault

Attachments:

bedload_wi...ecking.f (5KB)

The administrator has disabled public write access.

Wilcock and Crowe subroutine modification - segfault 1 year 6 months ago #42559

kopmann
OFFLINE
Senior Boarder
Posts: 106
Thank you received: 65

Dear André,

I tried your subroutine but avoid changing the parameterlist because then you have to change the calling subroutine as well.
With "use declarations_gaia" you can include global variables. Because dzfdx,dzfdy and unorm are unfortunately not global variables they have to be computed again in the subroutine. It is not the nicest/most efficient way to do but it works.
I tried the attached subroutine with the yen-example and the results look OK.

Best regards,
Rebekka

Attachments:

bedlad_wil...e_gaia.f (6KB)

The administrator has disabled public write access.

The following user(s) said Thank You: Renault

Wilcock and Crowe subroutine Recking modification 1 year 6 months ago #42560

Renault
OFFLINE
Senior Boarder
Posts: 120
Thank you received: 33

Hi Rebekka,

This is excellent, thank you so much! I was suspecting that I wasn't able to access the DZFD* and UNORM variables, but I didn't think of the problems in modifying the function call. Clearly, I still have a lot to learn with FORTRAN and GAIA... As much as there is to learn, I'm grateful it's possible to make these modifications as end users (with some help from knowledgeable staff

).

I will test this on my end and report back, although if it worked for you, it should work for me too.

Incidentally (for anyone else looking at this), this is one of the modifications to Wilcock and Crowe suggested in the GAIA user manual. I think I downloaded the Recking et al. paper based on the user manual's recommendation and much later decided to apply the correction, without thinking of the manual! Either way, a modified subroutine exists now