UNST-9784: Turn on 1d solver for 1D2D links again.

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_data/unstruc_model.f90

@@ -428,6 +428,7 @@ subroutine loadModel(filename)
       use m_realan, only: realan
       use m_filez, only: oldfil
       use unstruc_messages, only: threshold_abort
+      use m_default1d2d, only: set_default_1d2d_crs
       use m_readCrossSections, only: readCrossSectionDefinitions
 
       character(*), intent(inout) :: filename !< Name of file to be read (in current directory or with full path).
@@ -509,6 +510,7 @@ subroutine loadModel(filename)
       ! set administration arrays and fill cross section list. So getbobs for 1d can be called.
       call timstrt('Initialise 1d administration', timerHandle)
       call initialize_1dadmin(network, network%numl, numl)
+      call set_default_1d2d_crs(network)
       call timstop(timerHandle)
 
       if (getMaxErrorLevel() >= LEVEL_ERROR) then

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/compute/getprof_1d.f90

@@ -45,6 +45,7 @@ subroutine getprof_1D(L, hprL, area, width, japerim, calcConv, perim)
       use m_cross_helper, only: getconveyance, getcrossdischarge
       use m_flowtimes, only: time1, times_update_roughness
       use m_get_chezy, only: get_chezy
+      use m_roughness, only: getchezy
 
       integer :: L, japerim, calcConv
       real(kind=dp) :: hprL !< hoogte in profiel
@@ -122,7 +123,7 @@ subroutine getprof_1D(L, hprL, area, width, japerim, calcConv, perim)
 
          return
 
-      else if (abs(kcu(ll)) == 1 .and. network%loaded) then !flow1d used only for 1d channels and not for 1d2d roofs and gullies
+      else if ((abs(kcu(ll)) == 1 .or. abs(kcu(ll)) == 5) .and. network%loaded) then !flow1d used only for 1d channels and not for 1d2d roofs and gullies
          cz = 0.0_dp
 
          if (japerim == 0) then ! calculate total area and volume
@@ -134,17 +135,25 @@ subroutine getprof_1D(L, hprL, area, width, japerim, calcConv, perim)
             if (calcConv == 1) then
                u1L = u1(LL)
                q1L = q1(LL)
-               k1 = ln(1, LL)
-               k2 = ln(2, LL)
-               s1L = acl(L) * s1(k1) + (1.0_dp - acl(L)) * s1(k2)
                dpt = hu(L)
                cz = 0.0_dp
                if (network%rgs%timeseries_defined) then
                   factor = max(0.0_dp, (time1 - times_update_roughness(1)) / (times_update_roughness(2) - times_update_roughness(1)))
                else
                   factor = 1.0_dp
                end if
-               call getconveyance(network, dpt, u1L, q1L, s1L, LL, perim_sub, af_sub, conv, cz_sub, cz, area, perim, factor)
+               if (abs(kcu(ll)) == 1) then
+                  k1 = ln(1, LL)
+                  k2 = ln(2, LL)
+                  s1L = acl(L) * s1(k1) + (1.0_dp - acl(L)) * s1(k2)
+                  call getconveyance(network, dpt, u1L, q1L, s1L, LL, perim_sub, af_sub, conv, cz_sub, cz, area, perim, factor)
+               else ! 1D2Dlink
+                  frcn = network%crs%cross(network%adm%line2cross(LL, 2)%c1)%frictionvaluepos(1) !>
+                  friction_type = network%crs%cross(network%adm%line2cross(LL, 2)%c1)%frictiontypepos(1)
+                  cz = getchezy(friction_type, frcn, area / perim, dpt, u1(LL))
+                  cz_sub(1) = cz
+                  conv = cz_sub(1) * af_sub(1) * sqrt(af_sub(1) / perim_sub(1))
+               end if
 
                ! For sediment transport the discharge in the main channel is required:
                ! Qmain/ QT = Kmain/KT -> u_main = Kmain/KT * (AT/Amain)

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/compute/m_longculverts.f90

@@ -462,7 +462,6 @@ subroutine loadLongCulvertsAsNetwork(structurefile, jaKeepExisting, ierr)
       type(tree_data), pointer :: str_ptr
       character(len=IdLen) :: typestr
       character(len=IdLen) :: st_id
-      character(len=IdLen) :: csDefId
       character(len=IdLen) :: txt
       integer :: readerr, nstr, i, numcoords
       integer, allocatable, dimension(:) :: links
@@ -564,7 +563,7 @@ subroutine loadLongCulvertsAsNetwork(structurefile, jaKeepExisting, ierr)
                call prop_get(str_ptr, '', 'contactId', longculverts(nlongculverts)%contactID, success)
             end if
             if (success) then
-               call prop_get(str_ptr, '', 'csDefId', csDefId, success)
+               call prop_get(str_ptr, '', 'csDefId', longculverts(nlongculverts)%csDefId, success)
                if (.not. success) then
                   call SetMessage(LEVEL_ERROR, 'csDefId not found for long culvert: '//trim(st_id))
                end if
@@ -579,7 +578,7 @@ subroutine loadLongCulvertsAsNetwork(structurefile, jaKeepExisting, ierr)
             longculverts(nlongculverts)%netlinks = -999
 
             if (newculverts) then
-               call addlongculvertcrosssections(network, longculverts(nlongculverts)%branchid, csDefId, longculverts(nlongculverts)%bl, iref)
+               call addlongculvertcrosssections(network, longculverts(nlongculverts)%branchid, longculverts(nlongculverts)%csDefId, longculverts(nlongculverts)%bl, iref)
                if (iref > 0) then
                   ! Use top (#2) of tabulated cross section definition to derive width and height
                   longculverts(nlongculverts)%width = network%CSDefinitions%Cs(iref)%totalwidth(2)
@@ -702,6 +701,7 @@ end subroutine reallocLongCulverts
    subroutine longculvertsToProfs(skiplinks)
       use network_data
       use m_flowgeom
+      use unstruc_channel_flow, only: network
 
       logical, intent(in) :: skiplinks !< Skip determining the flow links or not
 
@@ -773,10 +773,8 @@ subroutine longculvertsToProfs(skiplinks)
             end do
             Lf = abs(longculverts(ilongc)%flowlinks(1))
             if (Lf > 0) then
+               call add_longculvert_1D2D_crosssection(network, Lf, longculverts(ilongc)%branchid, longculverts(ilongc)%csdefId)
                wu(Lf) = longculverts(ilongc)%width
-               prof1D(1, Lf) = wu(Lf)
-               prof1D(2, Lf) = longculverts(ilongc)%height
-               prof1D(3, Lf) = -2
                bob(1, Lf) = longculverts(ilongc)%bl(1)
                bob(2, Lf) = bl(ln(2, Lf))
             end if
@@ -1239,6 +1237,41 @@ subroutine addlongculvertcrosssections(network, branchId, csdefId, zpl, iref)
       end if
 
    end subroutine addlongculvertcrosssections
+
+   !> add special 1D2D crossection for the longculvert and add it to the line2cross array
+   subroutine add_longculvert_1D2D_crosssection(network, flowlink, branchId, csdefId)
+      use precision, only: dp
+      use m_hash_search
+      use m_readCrossSections
+      use m_network
+      type(t_network), intent(inout) :: network !< Network structure
+      integer, intent(in) :: flowlink !< Flowlink number on which to place the cross section. Should be 1D2D link belonging to the long culvert
+      character(len=IdLen), intent(in) :: branchId !< Branch id on which to place the cross section
+      character(len=IdLen), intent(in) :: csdefId !< Id of cross section definition
+
+      integer :: iref, icrs
+      integer :: indx
+      type(t_CrossSection), pointer :: pCrs
+      character(len=5) :: kchar
+
+      indx = hashsearch(network%brs%hashlist, branchId)
+      iref = hashsearch(network%CSDefinitions%hashlist, csdefId)
+      if (indx > 0 .and. iref > 0) then
+         if (network%crs%count + 1 > network%crs%size) then
+            call realloc(network%crs)
+         end if
+         icrs = network%crs%count + 1
+         pCrs => network%crs%cross(icrs)
+         write (kchar, '(I0)') flowlink
+         pCrs%csid = trim(branchId)//'_1D2D_'//trim(kchar)
+         call finalizeCrs(network, pCrs, iref, icrs)
+         network%adm%line2cross(flowlink, :)%c1 = icrs
+         network%adm%line2cross(flowlink, :)%c2 = icrs
+         network%adm%line2cross(flowlink, :)%f = 1.0_dp
+         network%adm%line2cross(flowlink, :)%distance = 0.0_dp
+      end if
+
+   end subroutine add_longculvert_1D2D_crosssection
    !> Add new branch iformation to the network. Only add necessary information for long culverts (incomplete!)
    subroutine add_longculvert_branch(network, longculvert)
       use precision, only: dp

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/default1d2d.f90

@@ -0,0 +1,96 @@
+module m_default1d2d
+   implicit none
+
+contains
+
+   subroutine set_default_1d2d_crs(network)
+      use precision, only: dp
+      use m_network, only: t_network
+      use m_flowgeom, only: wu1Duni5, hh1Duni5
+      use m_physcoef, only: ifrctypuni, frcuni1d2d
+      use m_CrossSections, only: AddCrossSectionDefinition, realloc
+      use network_data, only: numl1d, kn
+      use m_GlobalParameters, only: t_chainage2cross
+      use m_readcrosssections, only: finalizeCrs
+
+      type(t_network), target, intent(inout) :: network
+
+      ! Local variables
+      integer :: idef, icrs, L
+      integer :: numLevels
+      real(kind=dp) :: level(2)
+      real(kind=dp) :: width(2)
+      real(kind=dp) :: plains(3)
+      real(kind=dp) :: crestLevel, baseLevel, flowArea, totalArea
+      logical :: closed
+      logical :: groundlayerUsed
+      real(kind=dp) :: groundlayer
+      character(len=:), allocatable :: id
+      integer, parameter :: MANNING = 1
+      type(t_chainage2cross), dimension(:, :), allocatable :: temp_line2cross
+
+      ! Create a default 1D2D rectangular cross-section definition
+      id = 'default_1d2d_rect'
+      closed = .true.
+      groundlayerUsed = .false.
+      groundlayer = 0.0_dp
+
+      ! Set up rectangular profile with wu1Duni5 (width) and hh1Duni5 (height)
+      numLevels = 2 ! closed rectangular section has 2 levels
+
+      level(1) = 0.0_dp ! Bottom level
+      level(2) = hh1Duni5 ! Top level (height)
+      width(1) = wu1Duni5 ! Width at bottom
+      width(2) = wu1Duni5 ! Width at top (same for rectangle)
+
+      plains = 0.0_dp !wu1Duni5
+      crestLevel = 0.0_dp
+      baseLevel = 0.0_dp
+      flowArea = 0.0_dp
+      totalArea = 0.0_dp
+
+      ! Add the cross-section definition to the network
+      idef = AddCrossSectionDefinition(network%CSDefinitions, id, numLevels, level, width, &
+                                       width, plains, crestLevel, baseLevel, flowArea, totalArea, &
+                                       closed, groundlayerUsed, groundlayer)
+
+      ! Set up friction section using default 1D2D friction values
+      network%CSDefinitions%CS(idef)%frictionSectionsCount = 1
+      allocate (network%CSDefinitions%CS(idef)%frictionSectionID(1))
+      allocate (network%CSDefinitions%CS(idef)%frictionSectionIndex(1))
+      allocate (network%CSDefinitions%CS(idef)%frictionType(1))
+      allocate (network%CSDefinitions%CS(idef)%frictionValue(1))
+
+      network%CSDefinitions%CS(idef)%frictionSectionID(1) = ''
+      network%CSDefinitions%CS(idef)%frictionSectionIndex(1) = 0
+
+      ! Use default friction values for 1D2D links
+      ! For lateral 1D2D links, use frcuni1d2d with ifrctypuni
+      network%CSDefinitions%CS(idef)%frictionType(1) = ifrctypuni
+      network%CSDefinitions%CS(idef)%frictionValue(1) = frcuni1d2d
+
+      ! Now create a cross-section instance
+      if (network%crs%count + 1 > network%crs%size) then
+         call realloc(network%crs)
+      end if
+      icrs = network%crs%count + 1
+      call finalizeCrs(network, network%crs%cross(icrs), idef, icrs)
+      ! Note: For street inlets and roof gutters, you would use:
+      ! frictionType = MANNING
+      ! frictionValue = frcunistreetinlet or frcuniroofgutterpipe
+      temp_line2cross = network%adm%line2cross
+      deallocate (network%adm%line2cross)
+      allocate (network%adm%line2cross(numl1d, 3))
+      network%adm%line2cross(1:size(temp_line2cross, 1), :) = temp_line2cross
+      do L = 1, numl1d
+         if (kn(3, L) == 5 .and. all(network%adm%line2cross(L, :)%c1 == -1)) then !> don't overwrite long culvert lin2cross
+            network%adm%line2cross(L, :)%c1 = icrs
+            network%adm%line2cross(L, :)%c2 = icrs
+            network%adm%line2cross(L, :)%f = 1.0_dp
+            network%adm%line2cross(L, :)%distance = 0.0_dp
+         end if
+      end do
+
+   end subroutine set_default_1d2d_crs
+
+end module m_default1d2d

src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt

@@ -25,12 +25,21 @@ f90twtest(test_pol_to_cellmask
     VISUAL_STUDIO_PATH engines_gpl/dflowfm/test
 )
 
+f90twtest(test_getprof_1d
+    CFILES test_getprof_1d.cpp
+    F90FILES test_getprof_1d.f90
+    F2HFILES ${CMAKE_CURRENT_SOURCE_DIR}/test_getprof_1d.f90
+    OUTDIR "${CMAKE_CURRENT_BINARY_DIR}"
+    LIBRARIES dflowfm_kernel test_helpers
+    VISUAL_STUDIO_PATH engines_gpl/dflowfm/test
+)
+
 f90twtest(test_longculverts
     CFILES test_longculverts.cpp
     F90FILES test_longculverts.f90
     F2HFILES ${CMAKE_CURRENT_SOURCE_DIR}/test_longculverts.f90
     OUTDIR "${CMAKE_CURRENT_BINARY_DIR}"
-    LIBRARIES dflowfm_kernel
+    LIBRARIES dflowfm_kernel test_helpers
     VISUAL_STUDIO_PATH engines_gpl/dflowfm/test
 )
 

src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/helpers/CMakeLists.txt

@@ -1,2 +1,2 @@
-add_library(test_helpers STATIC file_helpers.f90)
-target_link_libraries(test_helpers dflowfm_kernel)
+add_library(test_helpers STATIC file_helpers.f90 network_helpers.f90)
+target_link_libraries(test_helpers dflowfm_kernel)