diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90 index 59247766..71cd7904 100644 --- a/source/src/sfincs_continuity.f90 +++ b/source/src/sfincs_continuity.f90 @@ -54,6 +54,8 @@ subroutine compute_water_levels_regular(dt,t) ! integer :: nm integer :: isrc + integer :: idrn, jin, jout, nmin, nmout + real*4 :: dv, dva, zold, a, ain, aout ! integer :: iwm ! @@ -67,7 +69,7 @@ subroutine compute_water_levels_regular(dt,t) real*4 :: qnum real*4 :: qndm real*4 :: factime - real*4 :: dvol + real*4 :: dvol ! if (snapwave) then ! need to compute filtered water levels for snapwave ! @@ -85,24 +87,72 @@ subroutine compute_water_levels_regular(dt,t) ! First discharges (don't do this parallel, as it's probably not worth it) ! if (nsrcdrn > 0) then - ! - !$acc loop - do isrc = 1, nsrcdrn - ! + ! + ! Plain point sources: apply the volume change and floor at the bed + ! + !$acc loop seq + do isrc = 1, nsrc + ! nm = nmindsrc(isrc) - ! + ! if (crsgeo) then - ! - zs(nmindsrc(isrc)) = max(zs(nm) + qtsrc(isrc) * dt / cell_area_m2(nm), zb(nm)) - ! + a = cell_area_m2(nm) else - ! - zs(nmindsrc(isrc)) = max(zs(nm) + qtsrc(isrc) * dt / cell_area(z_flags_iref(nm)), zb(nm)) - ! + a = cell_area(z_flags_iref(nm)) endif - ! + ! + dv = qtsrc(isrc) * dt + if (dv < 0.0) then + srcdrn_vol_withdrawn = srcdrn_vol_withdrawn - dv + srcdrn_vol_shortfall = srcdrn_vol_shortfall + max(-((zs(nm) - zb(nm)) * a + dv), 0.0) + endif + ! + zs(nm) = max(zs(nm) + dv / a, zb(nm)) + ! enddo - ! + ! + ! Drainage structures: the donor cell is floored at the bed and the receiving cell + ! gets exactly the volume that was actually removed, so the pair conserves mass. + ! + !$acc loop seq + do idrn = 1, ndrn + ! + jin = nsrc + idrn * 2 - 1 + jout = nsrc + idrn * 2 + nmin = nmindsrc(jin) + nmout = nmindsrc(jout) + ! + if (nmin > 0 .and. nmout > 0) then + ! + if (crsgeo) then + ain = cell_area_m2(nmin) + aout = cell_area_m2(nmout) + else + ain = cell_area(z_flags_iref(nmin)) + aout = cell_area(z_flags_iref(nmout)) + endif + ! + dv = qtsrc(jout) * dt ! volume moved from intake to outfall this step (negative = reverse flow) + ! + if (dv >= 0.0) then + zold = zs(nmin) + zs(nmin) = max(zold - dv / ain, zb(nmin)) + dva = max((zold - zs(nmin)) * ain, 0.0) + zs(nmout) = zs(nmout) + dva / aout + else + zold = zs(nmout) + zs(nmout) = max(zold + dv / aout, zb(nmout)) + dva = max((zold - zs(nmout)) * aout, 0.0) + zs(nmin) = zs(nmin) + dva / ain + endif + ! + srcdrn_vol_withdrawn = srcdrn_vol_withdrawn + abs(dv) + srcdrn_vol_shortfall = srcdrn_vol_shortfall + (abs(dv) - dva) + ! + endif + ! + enddo + ! endif ! !$omp parallel & @@ -273,6 +323,8 @@ subroutine compute_water_levels_subgrid(dt,t) ! integer :: nm integer :: isrc + integer :: idrn, jin, jout, nmin, nmout + real*4 :: dva, zold ! integer :: iwm integer :: ind @@ -311,12 +363,51 @@ subroutine compute_water_levels_subgrid(dt,t) if (nsrcdrn > 0) then ! !$acc serial present( z_volume, nmindsrc, qtsrc ) - do isrc = 1, nsrcdrn + ! + ! Plain point sources: apply the volume change and floor at zero + ! + do isrc = 1, nsrc ! nm = nmindsrc(isrc) ! - if ((z_volume(nm) >= 0) .or. ((qtsrc(isrc)<0.0) .and. (z_volume(nm) >= 0))) then - z_volume(nm) = z_volume(nm) + qtsrc(isrc) * dt + dv = qtsrc(isrc) * dt + if (dv < 0.0) then + srcdrn_vol_withdrawn = srcdrn_vol_withdrawn - dv + srcdrn_vol_shortfall = srcdrn_vol_shortfall + (max(z_volume(nm) + dv, 0.0) - (z_volume(nm) + dv)) + endif + z_volume(nm) = max(z_volume(nm) + dv, 0.0) + ! + enddo + ! + ! Drainage structures: the donor cell is floored at zero and the receiving cell + ! gets exactly the volume that was actually removed, so the pair conserves mass. + ! + do idrn = 1, ndrn + ! + jin = nsrc + idrn * 2 - 1 + jout = nsrc + idrn * 2 + nmin = nmindsrc(jin) + nmout = nmindsrc(jout) + ! + if (nmin > 0 .and. nmout > 0) then + ! + dv = qtsrc(jout) * dt ! volume moved from intake to outfall this step (negative = reverse flow) + ! + if (dv >= 0.0) then + zold = z_volume(nmin) + z_volume(nmin) = max(zold - dv, 0.0) + dva = max(zold - z_volume(nmin), 0.0) + z_volume(nmout) = max(z_volume(nmout) + dva, 0.0) + else + zold = z_volume(nmout) + z_volume(nmout) = max(zold + dv, 0.0) + dva = max(zold - z_volume(nmout), 0.0) + z_volume(nmin) = max(z_volume(nmin) + dva, 0.0) + endif + ! + srcdrn_vol_withdrawn = srcdrn_vol_withdrawn + abs(dv) + srcdrn_vol_shortfall = srcdrn_vol_shortfall + (abs(dv) - dva) + ! endif ! enddo diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90 index 825f7af5..184c2ac9 100644 --- a/source/src/sfincs_data.f90 +++ b/source/src/sfincs_data.f90 @@ -89,6 +89,7 @@ module sfincs_data real*4 horton_kr_kd real*4 btrelax real*4 structure_relax + real*4 drainage_volfrac ! max fraction of the donor cell volume a drain may remove per time step (0-1) real*4 wiggle_factor real*4 wiggle_threshold real*4 uvlim @@ -802,6 +803,8 @@ module sfincs_data real*4, dimension(:), allocatable :: drainage_distance integer*1, dimension(:), allocatable :: drainage_status real*4, dimension(:), allocatable :: drainage_fraction_open + real*8 :: srcdrn_vol_withdrawn ! volume nominally withdrawn by drains/negative sources (m3) + real*8 :: srcdrn_vol_shortfall ! part of that not actually moved because the donor cell ran dry (m3) real*4, dimension(:), allocatable :: xsrc real*4, dimension(:), allocatable :: ysrc !!! diff --git a/source/src/sfincs_discharges.f90 b/source/src/sfincs_discharges.f90 index 2d3a7262..16487082 100644 --- a/source/src/sfincs_discharges.f90 +++ b/source/src/sfincs_discharges.f90 @@ -370,7 +370,7 @@ subroutine update_discharges(t, dt, tloop) ! if (ndrn > 0) then ! - !$acc serial, present( z_volume, zs, zb, nmindsrc, qtsrc, drainage_type, drainage_params ) + !$acc serial, present( z_volume, zs, zb, nmindsrc, qtsrc, drainage_type, drainage_params ) copyin( drainage_volfrac ) do idrn = 1, ndrn ! jin = nsrc + idrn * 2 - 1 @@ -619,32 +619,31 @@ subroutine update_discharges(t, dt, tloop) ! end select ! - ! Add some relaxation - ! structure_relax in seconds => gives ratio between new and old discharge (default 10s) - ! - qq = 1.0 / (structure_relax / dt) * qq + (1.0 - (1.0 / (structure_relax / dt))) * -qtsrc(jin) - ! - ! Limit discharge based on available volume in cell (regular or subgrid) - ! + ! Limit discharge to a fraction of the available volume in the donor cell (regular or subgrid). + ! With drainage_volfrac < 1 the cell volume decays geometrically instead of hitting zero, + ! so the limited discharge converges smoothly to the inflow instead of switching on and off. if (subgrid) then ! if (qq > 0.0) then - qq = min(qq, max(z_volume(nmin), 0.0) / dt) + qq = min(qq, drainage_volfrac * max(z_volume(nmin), 0.0) / dt) else - qq = max(qq, -max(z_volume(nmout), 0.0) / dt) + qq = max(qq, -drainage_volfrac * max(z_volume(nmout), 0.0) / dt) endif ! else ! if (qq > 0.0) then - qq = min(qq, max((zs(nmin) - zb(nmin)) * cell_area(z_flags_iref(nmin)), 0.0) / dt) + qq = min(qq, drainage_volfrac * max((zs(nmin) - zb(nmin)) * cell_area(z_flags_iref(nmin)), 0.0) / dt) else - qq = max(qq, -max((zs(nmout) - zb(nmout)) * cell_area(z_flags_iref(nmout)), 0.0) / dt) + qq = max(qq, -drainage_volfrac * max((zs(nmout) - zb(nmout)) * cell_area(z_flags_iref(nmout)), 0.0) / dt) endif ! endif - ! - qtsrc(jin) = -qq + ! + ! Add some relaxation (applied LAST, so the applied discharge is a smooth low-pass). + qq = 1.0 / (structure_relax / dt) * qq + (1.0 - (1.0 / (structure_relax / dt))) * -qtsrc(jin) + ! + qtsrc(jin) = -qq qtsrc(jout) = qq ! endif diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90 index a68e6fd1..808c263d 100644 --- a/source/src/sfincs_input.f90 +++ b/source/src/sfincs_input.f90 @@ -189,6 +189,8 @@ subroutine read_sfincs_input() call read_real_input(500,'btrelax',btrelax,3600.0) call read_logical_input(500,'wiggle_suppression', wiggle_suppression, .true.) call read_real_input(500,'structure_relax',structure_relax,10.0) + call read_real_input(500,'drainage_volfrac',drainage_volfrac,0.5) + drainage_volfrac = min(max(drainage_volfrac, 0.05), 1.0) call read_real_input(500,'wiggle_factor',wiggle_factor,0.1) call read_real_input(500,'wiggle_threshold',wiggle_threshold,0.1) call read_real_input(500, 'uvlim', uvlim, 10.0) diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90 index 8598771f..1442084b 100644 --- a/source/src/sfincs_lib.f90 +++ b/source/src/sfincs_lib.f90 @@ -281,6 +281,8 @@ function sfincs_initialize() result(ierr) dt = 1.0e-6 ! First time step very small min_dt = 1.0e-6 ! First time step very small dtavg = 0.0 ! average time step + srcdrn_vol_withdrawn = 0.0d0 + srcdrn_vol_shortfall = 0.0d0 maxdepth = 999.0 ! maximum depth over time step maxmaxdepth = 0.0 ! maximum depth over entire simulation nt = 0 ! number of time steps @@ -774,6 +776,12 @@ function sfincs_finalize() result(ierr) write(logstr,'(a,20f10.3)') ' Average time step (s) : ', dtavg call write_log(logstr, 1) ! + if (nsrcdrn > 0 .and. srcdrn_vol_withdrawn > 0.0d0) then + write(logstr,'(a,f14.1,a,f12.1,a,f7.3,a)') ' Drain/source withdrawal: ', srcdrn_vol_withdrawn, ' m3, not removed (dry cell): ', & + srcdrn_vol_shortfall, ' m3 (', 100.0d0 * srcdrn_vol_shortfall / srcdrn_vol_withdrawn, '%)' + call write_log(logstr, 1) + endif + ! call write_log('', 1) ! if (timestep_analysis) then