diff --git a/.gitignore b/.gitignore
index 04acf563a..55ffd2c50 100644
--- a/.gitignore
+++ b/.gitignore
@@ -65,6 +65,7 @@ source/third_party_open/netcdf/x64
source/sfincs/sfincs.opt.yaml
/source/sfincs_lib/*.yaml
/source/third_party_open/netcdf/netcdf-fortran-4.6.1/Debug
+source/build.log
/source/build_nvfortran_gpu_h7.sh
/source/build_nvfortran_gpu.sh
/source/Singularityfile-gpu.def
diff --git a/source/sfincs_lib/sfincs_lib.vfproj b/source/sfincs_lib/sfincs_lib.vfproj
index 9e9cac55a..cb8d1aa11 100644
--- a/source/sfincs_lib/sfincs_lib.vfproj
+++ b/source/sfincs_lib/sfincs_lib.vfproj
@@ -30,16 +30,12 @@
-
-
-
-
+
+
-
-
-
-
+
+
@@ -57,8 +53,7 @@
-
-
+
@@ -99,19 +94,15 @@
-
-
-
+
-
-
-
+
@@ -121,6 +112,7 @@
+
@@ -159,8 +151,7 @@
-
-
+
@@ -172,8 +163,7 @@
-
-
+
diff --git a/source/src/Makefile.am b/source/src/Makefile.am
index 1aa72660e..8b8bd1d79 100644
--- a/source/src/Makefile.am
+++ b/source/src/Makefile.am
@@ -88,12 +88,12 @@ libsfincs_la_SOURCES = \
sfincs_snapwave.f90 \
../third_party_open/utils/deg2utm.f90 \
sfincs_meteo.f90 \
- ../third_party_open/bicgstab/bicgstab_solver_ilu.f90 \
sfincs_nonhydrostatic.f90 \
sfincs_ncoutput_helpers.F90 \
sfincs_ncoutput.F90 \
sfincs_output.f90 \
sfincs_momentum.f90 \
+ sfincs_momentum_velocity.f90 \
sfincs_wavemaker.f90 \
sfincs_bathtub.f90 \
sfincs_openacc.f90 \
diff --git a/source/src/sfincs.f90 b/source/src/sfincs.f90
index 38fe30bb2..bab831811 100644
--- a/source/src/sfincs.f90
+++ b/source/src/sfincs.f90
@@ -15,8 +15,9 @@ program sfincs
!
! Set BMI flags to false
!
- bmi = .false.
+ bmi = .false.
use_qext = .false.
+ use_dzbext = .false.
!
ierr = sfincs_initialize()
!
diff --git a/source/src/sfincs_bmi.f90 b/source/src/sfincs_bmi.f90
index 388bcefa2..925aa1c38 100644
--- a/source/src/sfincs_bmi.f90
+++ b/source/src/sfincs_bmi.f90
@@ -23,7 +23,7 @@ module sfincs_bmi
public :: get_end_time
public :: get_time_step
public :: get_current_time
- public :: update_zbuv
+ public :: bmi_update_bed_level
public :: update_apparent_roughness
public :: get_sfincs_cell_index
public :: get_sfincs_cell_indices
@@ -119,6 +119,12 @@ function get_value_ptr(c_var_name, c_data) result(ierr) &
c_data = c_loc(subgrid_z_zmin)
case("qext")
c_data = c_loc(qext)
+ case("dzbext")
+ if (.not. allocated(dzbext)) then
+ allocate(dzbext(np))
+ dzbext = 0.0
+ endif
+ c_data = c_loc(dzbext)
case("uorb")
c_data = c_loc(uorb)
case default
@@ -150,6 +156,12 @@ function get_var_shape(c_var_name, var_shape) result(ierr) &
var_shape(1) = size(z_index_z_n)
case("qext")
var_shape(1) = size(qext)
+ case("dzbext")
+ if (.not. allocated(dzbext)) then
+ allocate(dzbext(np))
+ dzbext = 0.0
+ endif
+ var_shape(1) = size(dzbext)
case default
ierr = -1
end select
@@ -170,7 +182,7 @@ function get_var_type(c_var_name, c_type) result(ierr) &
var_name = char_array_to_string(c_var_name, strlen(c_var_name, BMI_LENVARADDRESS))
select case(var_name)
- case("z_xz", "z_yz", "zb", "subgrid_z_zmin", "qext", "uorb")
+ case("z_xz", "z_yz", "zb", "subgrid_z_zmin", "qext", "dzbext", "uorb")
type_name = "float"
case("zs")
type_name = "double"
@@ -198,7 +210,7 @@ function get_var_rank(c_var_name, rank) result(ierr) &
var_name = char_array_to_string(c_var_name, strlen(c_var_name, BMI_LENVARADDRESS))
select case(var_name)
- case("z_xz", "z_yz", "zs", "zb", "subgrid_z_zmin", "qext", "uorb")
+ case("z_xz", "z_yz", "zs", "zb", "subgrid_z_zmin", "qext", "dzbext", "uorb")
rank = 1
case default
ierr = -1
@@ -227,7 +239,23 @@ function set_logical(c_flag_name, ival) result(ierr) bind(C, name="set_logical")
select case(flag_name)
case("qext")
use_qext = bval
- !write(*,*)'use_qext = ', use_qext
+ !write(*,*)'use_qext = ', use_qext
+ case("dzbext")
+ !
+ ! Lazily allocate the external delta-bed-level array on first enable.
+ ! Once allocated we keep it around; toggling the flag off later just
+ ! disables the update path without freeing memory (same pattern as
+ ! qext is handled elsewhere).
+ !
+ if (bval .and. .not. allocated(dzbext)) then
+ !
+ allocate(dzbext(np))
+ dzbext = 0.0
+ !
+ endif
+ !
+ use_dzbext = bval
+ !
case default
ierr = -1
end select
@@ -274,14 +302,15 @@ function get_current_time(tcurrent) result(ierr) bind(C, name="get_current_time"
end function get_current_time
- function update_zbuv() result(ierr) bind(C, name="update_zbuv")
- ! Update bed level at uv points
- !DEC$ ATTRIBUTES DLLEXPORT :: update_zbuv
+ function bmi_update_bed_level() result(ierr) bind(C, name="update_bed_level")
+ ! Apply dzbext to the bed-level arrays (zb or subgrid_z_zmin/zmax and
+ ! subgrid_uv_zmin/zmax) and rebuild zbuvmx for non-subgrid runs.
+ !DEC$ ATTRIBUTES DLLEXPORT :: bmi_update_bed_level
integer(kind=c_int) :: ierr
- call compute_zbuvmx()
+ call update_bed_level()
ierr = 0
-
- end function update_zbuv
+
+ end function bmi_update_bed_level
function update_apparent_roughness() result(ierr) bind(C, name="update_apparent_roughness")
! Update apparent roughness at uv points
diff --git a/source/src/sfincs_continuity.f90 b/source/src/sfincs_continuity.f90
index 592477666..24490089c 100644
--- a/source/src/sfincs_continuity.f90
+++ b/source/src/sfincs_continuity.f90
@@ -569,6 +569,15 @@ subroutine compute_water_levels_subgrid(dt,t)
!
endif
!
+ ! Effective bed for subgrid models with velocity-form advection and/or
+ ! non-hydrostatics: zb = zs - V/A, the bed consistent with the volume
+ ! continuity conserves (the file zb is not a valid conveyance bed).
+ ! Bed SLOPES elsewhere stay frozen at their initialization values.
+ !
+ if (zb_effective) then
+ zb(nm) = zs(nm) - max(z_volume(nm) / a, 0.0)
+ endif
+ !
!
if (wiggle_suppression) then
!
@@ -618,9 +627,9 @@ subroutine compute_water_levels_subgrid(dt,t)
enddo
!$omp end do
!$omp end parallel
- !
+ !
!$acc end parallel
- !
+ !
end subroutine
subroutine compute_store_variables(dt)
diff --git a/source/src/sfincs_data.f90 b/source/src/sfincs_data.f90
index 534757471..a409ca979 100644
--- a/source/src/sfincs_data.f90
+++ b/source/src/sfincs_data.f90
@@ -16,6 +16,7 @@ module sfincs_data
!!!
logical :: bmi
logical :: use_qext
+ logical :: use_dzbext
!!!
!!! Constants
!!!
@@ -97,10 +98,7 @@ module sfincs_data
!real*4 dzdsbnd
!real*4 manningbnd
real*4 nuviscfac ! Factor on viscosity for 'difficult' points. Used in sfincs_momentum.f90.
- real*4 nh_fnudge
- real*4 nh_tstop
- integer nh_itermax
- real*4 nh_tol
+ logical zb_effective ! subgrid + (velocity scheme or nonh): zb is recomputed every step in continuity as the effective bed zs - z_volume/area
real*4 runup_gauge_depth
real*4 factor_wind
real*4 factor_pres
@@ -114,6 +112,7 @@ module sfincs_data
integer cd_nr
integer baro
integer advection_scheme
+ integer momentum_scheme ! 0 = Bates flux form (default), 1 = velocity form
!
character*256 :: depfile
character*256 :: mskfile
@@ -336,7 +335,6 @@ module sfincs_data
integer*1, dimension(:), allocatable :: kfuv
integer*1, dimension(:), allocatable :: mask_adv
integer*1, dimension(:), allocatable :: scs_rain ! logic if previous time step was raining
- integer*1, dimension(:), allocatable :: mask_nonh
!
! Quadtree
!
@@ -597,6 +595,7 @@ module sfincs_data
real*4, dimension(:), allocatable :: zs0
real*4, dimension(:), allocatable :: zsderv
real*4, dimension(:), allocatable, target :: qext
+ real*4, dimension(:), allocatable, target :: dzbext
real*4, dimension(:), allocatable, target :: uorb
real*4, dimension(:), allocatable :: gnapp2
!
@@ -1018,6 +1017,7 @@ subroutine finalize_parameters()
if(allocated(uv0)) deallocate(uv0)
if(allocated(twet)) deallocate(twet)
if(allocated(qext)) deallocate(qext)
+ if(allocated(dzbext)) deallocate(dzbext)
!
! if(allocated(huu)) deallocate(huu)
! if(allocated(hvv)) deallocate(hvv)
diff --git a/source/src/sfincs_domain.f90 b/source/src/sfincs_domain.f90
index c06595d59..4948b2268 100644
--- a/source/src/sfincs_domain.f90
+++ b/source/src/sfincs_domain.f90
@@ -145,6 +145,7 @@ subroutine initialize_mesh()
! Initialize SFINCS domain (indices, flags and neighbors)
!
use sfincs_data
+ use sfincs_nonhydrostatic, only: mask_nonh
use quadtree
!
implicit none
@@ -1274,7 +1275,12 @@ subroutine initialize_mesh()
!
if (kcs(nm) /= 1) cycle ! not a regular point
!
- ! Check if point has 4 neighbors with kcs = 1
+ ! Check that the point has regular (kcs = 1) neighbors. In the general
+ ! 2D case all four (left/right/bottom/top) are required. For a true 1-D
+ ! model (nmax == 1) cells legitimately have no bottom/top neighbor, so
+ ! only left/right are required there; the non-hydrostatic solver then
+ ! degrades to a 1-D (x-only) pressure Poisson problem, which all the
+ ! downstream assembly already supports through its per-direction guards.
!
! Left
!
@@ -1284,7 +1290,7 @@ subroutine initialize_mesh()
!
if (nmd == 0) cycle ! no neighbor
!
- if (kcs(nmd) /= 1) cycle ! neighbor is not a regular point
+ if (kcs(nmd) < 1 .or. kcs(nmd) > 3) cycle ! allow regular (1) or open-boundary (2/3) neighbour
!
! Right
!
@@ -1294,29 +1300,33 @@ subroutine initialize_mesh()
!
if (nmu == 0) cycle ! no neighbor
!
- if (kcs(nmu) /= 1) cycle ! neighbor is not a regular point
+ if (kcs(nmu) < 1 .or. kcs(nmu) > 3) cycle ! allow regular (1) or open-boundary (2/3) neighbour
!
- ! Bottom
- !
- if (z_flags_nd(nm) /= 0) cycle ! not a regular neighbor
- !
- ndm = z_index_z_nd1(nm)
- !
- if (ndm == 0) cycle ! no neighbor
- !
- if (kcs(ndm) /= 1) cycle ! neighbor is not a regular point
- !
- ! Top
- !
- if (z_flags_nu(nm) /= 0) cycle ! not a regular neighbor
- !
- num = z_index_z_nu1(nm)
- !
- if (num == 0) cycle ! no neighbor
- !
- if (kcs(num) /= 1) cycle ! neighbor is not a regular point
+ if (nmax > 1) then
+ !
+ ! Bottom
+ !
+ if (z_flags_nd(nm) /= 0) cycle ! not a regular neighbor
+ !
+ ndm = z_index_z_nd1(nm)
+ !
+ if (ndm == 0) cycle ! no neighbor
+ !
+ if (kcs(ndm) < 1 .or. kcs(ndm) > 3) cycle ! allow regular (1) or open-boundary (2/3) neighbour
+ !
+ ! Top
+ !
+ if (z_flags_nu(nm) /= 0) cycle ! not a regular neighbor
+ !
+ num = z_index_z_nu1(nm)
+ !
+ if (num == 0) cycle ! no neighbor
+ !
+ if (kcs(num) < 1 .or. kcs(num) > 3) cycle ! allow regular (1) or open-boundary (2/3) neighbour
+ !
+ endif
!
- ! This cell has 4 regular neighbors, so copy from quadtree mask
+ ! This cell has the required regular neighbors, so copy from quadtree mask
!
mask_nonh(nm) = quadtree_nonh_mask(index_quadtree_in_sfincs(nm))
!
@@ -1715,25 +1725,27 @@ subroutine initialize_bathymetry()
!
call read_subgrid_file()
!
- ! In case of nonh, we also need zb
+ ! In case of nonh or the velocity-form momentum scheme, we also need a per-cell
+ ! bed level. Subgrid models have no reliable file zb, so use the EFFECTIVE bed
+ ! implied by the subgrid tables at full wetness, zb = z_zmax - z_volmax/area:
+ ! the cell-mean bed. From the first time step onward continuity maintains the
+ ! water-level-dependent effective bed zb = zs - z_volume/area (zb_effective).
+ ! The non-hydrostatic initialization freezes its w_b bed slopes from the
+ ! cell-mean bed computed here.
!
- if (nonhydrostatic) then
+ if (nonhydrostatic .or. momentum_scheme == 1) then
!
allocate(zb(np))
!
- if (use_quadtree) then
- !
- do ip = 1, np
- zb(ip) = quadtree_zz(index_quadtree_in_sfincs(ip))
- enddo
- !
- else
+ do ip = 1, np
!
- ! Produce error message
- !
- call write_log('Error! : combination of nonhydrostatic solver with quadtree grid with nr_levels > 1 is not supported ', 1)
+ if (crsgeo) then
+ zb(ip) = subgrid_z_zmax(ip) - subgrid_z_volmax(ip) / cell_area_m2(ip)
+ else
+ zb(ip) = subgrid_z_zmax(ip) - subgrid_z_volmax(ip) / cell_area(z_flags_iref(ip))
+ endif
!
- endif
+ enddo
!
endif
!
@@ -1743,24 +1755,100 @@ subroutine initialize_bathymetry()
!
end subroutine
- subroutine compute_zbuvmx()
+ subroutine update_bed_level()
+ !
+ ! Apply the externally-supplied delta bed level (dzbext) to the kernel's
+ ! bed-level arrays and refresh derived quantities at uv points.
+ !
+ ! Non-subgrid mode:
+ ! zb = zb + dzbext (cell centres)
+ ! zbuvmx(ip) = max(zb(nm), zb(nmu)) + huthresh (uv points, rebuilt)
+ !
+ ! Subgrid mode:
+ ! subgrid_z_zmin/zmax shift rigidly by dzbext at the cell centre.
+ ! subgrid_uv_zmin/zmax shift rigidly by the average dzbext of the two
+ ! neighbouring cells, with subgrid_uv_zmin clamped from below by the
+ ! larger of the two updated cell-centre subgrid_z_zmin values so the uv
+ ! minimum can never sit below either neighbour's minimum.
+ !
+ ! The caller (Python via BMI) owns the lifecycle of dzbext: this routine
+ ! does not zero it out after applying. When use_dzbext is .false. the
+ ! routine still rebuilds zbuvmx in non-subgrid mode (cheap, and matches the
+ ! historical behaviour of compute_zbuvmx).
!
use sfincs_data
!
integer :: ip
integer :: nm
integer :: nmu
+ integer :: ilevel
+ real*4 :: avg_dzb
!
- do ip = 1, npuv
+ if (.not. subgrid) then
!
- nm = uv_index_z_nm(ip)
- nmu = uv_index_z_nmu(ip)
+ ! Non-subgrid path: shift zb, then rebuild zbuvmx for every uv point.
!
- zbuvmx(ip) = max(zb(nm), zb(nmu)) + huthresh
- !
- enddo
+ if (use_dzbext) then
+ !
+ zb(:) = zb(:) + dzbext(:)
+ !
+ endif
+ !
+ do ip = 1, npuv
+ !
+ nm = uv_index_z_nm(ip)
+ nmu = uv_index_z_nmu(ip)
+ !
+ zbuvmx(ip) = max(zb(nm), zb(nmu)) + huthresh
+ !
+ enddo
+ !
+ else
+ !
+ ! Subgrid path: only do anything when an external delta has been set.
+ !
+ if (use_dzbext) then
+ !
+ ! Cell-centre arrays shift rigidly with dzbext.
+ !
+ subgrid_z_zmin(:) = subgrid_z_zmin(:) + dzbext(:)
+ subgrid_z_zmax(:) = subgrid_z_zmax(:) + dzbext(:)
+ !
+ ! The per-level volume->water-level table holds ABSOLUTE levels, so it
+ ! must translate with the bed as well -- otherwise partially-wet cells
+ ! (0 < z_volume < volmax, the interpolation branch in
+ ! compute_water_levels_subgrid) would not follow the moving bed. A rigid
+ ! vertical shift leaves the volume bins (volmax / dzvol) unchanged, so
+ ! only the levels move. Keep the -20 m floor used in the level lookups.
+ !
+ do ilevel = 1, subgrid_nlevels
+ subgrid_z_dep(ilevel, :) = max(subgrid_z_dep(ilevel, :) + dzbext(:), -20.0)
+ enddo
+ !
+ ! UV-point arrays shift by the average delta of the two neighbours,
+ ! then clamp uv_zmin from below by the higher of the two updated
+ ! cell-centre minima.
+ !
+ do ip = 1, npuv
+ !
+ nm = uv_index_z_nm(ip)
+ nmu = uv_index_z_nmu(ip)
+ !
+ avg_dzb = 0.5 * (dzbext(nm) + dzbext(nmu))
+ !
+ subgrid_uv_zmin(ip) = subgrid_uv_zmin(ip) + avg_dzb
+ subgrid_uv_zmax(ip) = subgrid_uv_zmax(ip) + avg_dzb
+ !
+ subgrid_uv_zmin(ip) = max(subgrid_uv_zmin(ip), &
+ max(subgrid_z_zmin(nm), subgrid_z_zmin(nmu)))
+ !
+ enddo
+ !
+ endif
+ !
+ endif
!
- end subroutine
+ end subroutine update_bed_level
subroutine initialize_boundaries()
!
@@ -2309,7 +2397,10 @@ subroutine initialize_hydro()
t_zsmax = -999.0
endif
!
- uv = 0.0
+ ! NOTE: do NOT reset uv here. set_initial_conditions() (called above) computes the initial
+ ! velocity uv = q/hu from the restart/initial flux. The flux-form (Bates) scheme carries q and
+ ! recomputes uv every step, so this stray reset was harmless there, but the velocity-form
+ ! (momentum_scheme=velocity) scheme carries uv0 and needs the initial velocity preserved.
!
if (wind) then
!
diff --git a/source/src/sfincs_input.f90 b/source/src/sfincs_input.f90
index 40dc39105..bf81028e9 100644
--- a/source/src/sfincs_input.f90
+++ b/source/src/sfincs_input.f90
@@ -7,6 +7,7 @@ subroutine read_sfincs_input()
! Reads sfincs.inp
!
use sfincs_data
+ use sfincs_nonhydrostatic ! nonh_* input parameters live in the solver module
use sfincs_date
use sfincs_log
use sfincs_error
@@ -41,8 +42,9 @@ subroutine read_sfincs_input()
logical ok
!
character*256 wmsigstr
- character*256 advstr
- !
+ character*256 advstr
+ character*256 momstr
+ !
ok = check_file_exists('sfincs.inp', 'SFINCS input file', .true.)
!
open(500, file='sfincs.inp')
@@ -175,7 +177,8 @@ subroutine read_sfincs_input()
call read_logical_input(500, 'wavemaker_hinc', wavemaker_hinc, .false.) ! wavemaker include incident waves
!
! Numerical parameters
- call read_char_input(500,'advection_scheme',advstr,'upw1')
+ call read_char_input(500,'advection_scheme',advstr,'upw1')
+ call read_char_input(500,'momentum_scheme',momstr,'bates')
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)
@@ -188,12 +191,25 @@ subroutine read_sfincs_input()
! call read_real_input(500, 'dzdsbnd', dzdsbnd, 0.0001)
! call read_real_input(500, 'manningbnd', manningbnd, 0.024)
call read_real_input(500, 'nuviscfac', nuviscfac, 100.0)
- call read_logical_input(500, 'nonh', nonhydrostatic, .false.)
- call read_real_input(500, 'nh_fnudge', nh_fnudge, 0.9)
- call read_real_input(500, 'nh_tstop', nh_tstop, -999.0)
- call read_real_input(500, 'nh_tol', nh_tol, 0.001)
- call read_int_input(500, 'nh_itermax', nh_itermax, 100)
- call read_logical_input(500, 'h73table', h73table, .false.)
+ call read_logical_input(500, 'nonhydrostatic', nonhydrostatic, .false.)
+ call read_real_input(500, 'nonh_fnudge', nonh_fnudge, 1.0)
+ call read_real_input(500, 'nonh_tstop', nonh_tstop, -999.0)
+ call read_real_input(500, 'nonh_tol', nonh_tol, 0.01)
+ call read_int_input(500, 'nonh_itermax', nonh_itermax, 100)
+ call read_logical_input(500, 'nonh_movingbed', nonh_movingbed, .false.) ! add d(zb)/dt to the bottom kinematic w_b (moving-seafloor source via dzbext)
+ call read_real_input(500, 'nonh_filter', nonh_filter, 0.5) ! spatial 2dx filter on pnh (0 = off, ~0.25-0.5 damps grid mode)
+ call read_real_input(500, 'nonh_dzbmax', nonh_dzbmax, 0.1) ! cap on |d(zb)/dx| in bottom kinematic wb (default 0.1; clips near-vertical walls, leaves real slopes); 0 = no cap
+ call read_int_input(500, 'nonh_fadein', nonh_fadein, 0) ! open-boundary nonh fade-in width (cells): nonh ramps 0->full over N cells from the boundary; 0 = off
+ call read_real_input(500, 'nonh_brsteep', nonh_brsteep, 0.4) ! HFA breaking onset: nonh starts reducing when dzdt (=-d(hu)/dx) > nonh_brsteep*sqrt(g*h); 0 = off (XBeach default 0.4)
+ call read_real_input(500, 'nonh_brfr', nonh_brfr, 0.5) ! OPTIONAL NEOWAVE Froude breaking criterion: pnh=0 when |U|/sqrt(g*D) > nonh_brfr (~0.5), release < 0.3*nonh_brfr (~0.15); 0 = off -> use nonh_brsteep instead
+ call read_int_input(500, 'nonh_brsmooth', nonh_brsmooth, 0) ! breaking-flag smoothing passes ([1 2 1]/4 over nonh neighbours): ramps pnh out over ~brsmooth+1 cells at the breaking-zone edges instead of one face; 0 = sharp
+ call read_int_input(500, 'nonh_slsmooth', nonh_slsmooth, 0) ! frozen-bed-slope smoothing passes ([1 2 1]/4): bounds bed curvature d2zb/dx2 at slope breaks (island toe) -> suppresses trailing waves; 0 = off
+ call read_real_input(500, 'nonh_treform', nonh_treform, 1.0) ! breaking reformation time scale (s): released cells recover the nonh pressure gradually (brfac += dt/treform); 0 = instant recovery
+ call read_real_input(500, 'nonh_smoothbnd', nonh_smoothbnd, 0.5) ! strength of localized 2dx pnh smoothing in the fade-in zone (weight at boundary, ramps to 0 over the fade-in) and shallow zone; 0 = off
+ call read_real_input(500, 'nonh_smoothdep', nonh_smoothdep, 0.0) ! depth (m) below which the localized pnh smoothing also acts (shallow run-up / wall 2dx noise); weight ramps from full at D=0 to 0 at D=nonh_smoothdep; 0 = off
+ call read_real_input(500, 'nonh_disp', nonh_disp, 1.0) ! Keller-box vertical factor (default 1.0 = best dispersion, c(k) flat to Airy ~kd 2.5); 2.0 = strict linear-pressure single layer
+ call read_real_input(500, 'nonh_pmax', nonh_pmax, 0.0) ! depth limiter: cap |pnh| <= nonh_pmax*rho*g*H post-solve (caps wall/breaking/2dx spikes, leaves resolved waves untouched); 0 = off, ~1-2 typical
+ call read_logical_input(500, 'h73table', h73table, .false.)
call read_real_input(500, 'rugdepth', runup_gauge_depth, 0.05)
call read_logical_input(500, 'wave_enhanced_roughness', wave_enhanced_roughness, .false.)
call read_logical_input(500, 'use_bcafile', use_bcafile, .true.)
@@ -665,14 +681,41 @@ subroutine read_sfincs_input()
endif
endif
!
+ ! Momentum scheme : Bates flux form (default) or form (default)
+ !
+ if (trim(momstr) == 'bates') then
+ momentum_scheme = 0
+ call write_log('Info : momentum scheme : Bates (flux form)', 0)
+ else
+ momentum_scheme = 1
+ if (trim(momstr) /= 'velocity') then
+ write(logstr,*)'Warning : momentum scheme ', trim(momstr), ' not recognized! Using default velocity instead!'
+ call write_log(logstr, 1)
+ else
+ call write_log('Info : momentum scheme : velocity form', 0)
+ endif
+ endif
+ !
+ ! Effective bed level. Subgrid models have no reliable per-cell zb (the file
+ ! value is just the quadtree cell elevation), but the velocity-form advection
+ ! and the non-hydrostatic solver need one. For those combinations zb is
+ ! recomputed every step in sfincs_continuity as zs - z_volume/area: the bed
+ ! consistent with the volume continuity actually conserves. Bed SLOPES stay
+ ! frozen at their initialization (file zb) values (see sfincs_nonhydrostatic).
+ !
+ zb_effective = subgrid .and. (momentum_scheme == 1 .or. nonhydrostatic)
+ if (zb_effective) then
+ call write_log('Info : subgrid: using effective bed level zb = zs - V/A (velocity scheme / nonh)', 0)
+ endif
+ !
if (advection) then
!
! Make 1st order upwind the default scheme
- !
+ !
advection_scheme = 1
!
call write_log('Info : turning on advection', 0)
- !
+ !
if (trim(advstr) == 'original') then
advection_scheme = 0
call write_log('Info : advection scheme : Original', 0)
@@ -688,17 +731,17 @@ subroutine read_sfincs_input()
!
if (nonhydrostatic) then
!
- if (nh_tstop > 0.0) then
+ if (nonh_tstop > 0.0) then
!
! tstopnonh is provided so set it with respect to model reference time
!
- nh_tstop = t0 + nh_tstop
+ nonh_tstop = t0 + nonh_tstop
!
else
!
! tstopnonh is not provided so set it to tstop time + 999.0 s
!
- nh_tstop = t1 + 999.0
+ nonh_tstop = t1 + 999.0
!
endif
!
diff --git a/source/src/sfincs_lib.f90 b/source/src/sfincs_lib.f90
index a2fd05da9..c7267e019 100644
--- a/source/src/sfincs_lib.f90
+++ b/source/src/sfincs_lib.f90
@@ -19,6 +19,7 @@ module sfincs_lib
use sfincs_ncinput
use sfincs_ncoutput
use sfincs_momentum
+ use sfincs_momentum_velocity, only : compute_fluxes_velocity
use sfincs_continuity
use sfincs_snapwave
use sfincs_wavemaker
@@ -171,11 +172,11 @@ function sfincs_initialize() result(ierr)
!
! Initialize non-hydrostatic solver
!
- call write_log('Initialize non-hydrostatic solver ...', 0)
+ call write_log('Initialize non-hydrostatic solver ...', 0)
!
call initialize_nonhydrostatic()
!
- endif
+ endif
!
if (wavemaker) then
!
@@ -581,7 +582,19 @@ function sfincs_update(dtrange) result(ierr)
!
! First compute fluxes
!
- call compute_fluxes(dt, tloopflux)
+ if (momentum_scheme == 1) then
+ !
+ ! Momentum scheme 1 (velocity based)
+ !
+ call compute_fluxes_velocity(dt, tloopflux)
+ !
+ else
+ !
+ ! Original momentum scheme (flux based)
+ !
+ call compute_fluxes(dt, tloopflux)
+ !
+ endif
!
if (timestep_analysis) then
!
@@ -603,13 +616,13 @@ function sfincs_update(dtrange) result(ierr)
!
if (nonhydrostatic) then
!
- if (t < nh_tstop) then ! Check if non-hydrostatic corrections still need to be made
+ if (t < nonh_tstop) then ! Check if non-hydrostatic corrections still need to be made
!
! Apply non-hydrostatic pressure corrections to q and uv
!
call compute_nonhydrostatic(dt, tloopnonh)
!
- endif
+ endif
!
endif
!
@@ -746,6 +759,11 @@ function sfincs_finalize() result(ierr)
if (nonhydrostatic) then
write(logstr,'(a,f10.3,a,f5.1,a)') ' Time in non-hydrostatic: ', tloopnonh, ' (', 100 * tloopnonh / (tfinish_all - tstart_all), '%)'
call write_log(logstr, 1)
+ if (nh_solve_count > 0) then
+ write(logstr,'(a,f8.2,a,i0,a,i0,a)') ' CG iterations : ', &
+ real(nh_iter_total) / real(nh_solve_count), ' avg / ', nh_iter_max, ' max (', nh_solve_count, ' solves)'
+ call write_log(logstr, 1)
+ endif
endif
!
if (nrstructures>0) then
diff --git a/source/src/sfincs_momentum_velocity.f90 b/source/src/sfincs_momentum_velocity.f90
new file mode 100644
index 000000000..38300f97a
--- /dev/null
+++ b/source/src/sfincs_momentum_velocity.f90
@@ -0,0 +1,750 @@
+module sfincs_momentum_velocity
+ !
+ use sfincs_data
+ !
+ implicit none
+ !
+contains
+ !
+ subroutine compute_fluxes_velocity(dt, tloop)
+ !
+ ! Computes fluxes over subgrid u and v points
+ !
+ integer :: count0
+ integer :: count1
+ integer :: count_rate
+ integer :: count_max
+ real :: tloop
+ !
+ real*4 :: dt
+ !
+ integer :: ip
+ integer :: nm
+ integer :: nmu
+ !
+ integer :: idir
+ integer :: iref
+ integer :: itype
+ integer :: iuv
+ integer :: icuv
+ !
+ real*4 :: hu
+ real*4 :: dxuvinv
+ real*4 :: dxuv2inv
+ real*4 :: dyuvinv
+ real*4 :: dyuv2inv
+ real*4 :: adv
+ real*4 :: fcoriouv
+ real*4 :: frc
+ !
+ real*4 :: ufr
+ !
+ real*4 :: uu_nm
+ real*4 :: uu_nmd
+ real*4 :: uu_nmu
+ real*4 :: uu_ndm
+ real*4 :: uu_num
+ real*4 :: vu
+ !
+ real*4 :: zsu
+ real*4 :: dzuv
+ real*4 :: facint
+ real*4 :: gnavg2
+ real*4 :: fwmax
+ real*4 :: zmax
+ real*4 :: zmin
+ !
+ real*4 :: dqxudx
+ real*4 :: dqyudy
+ real*4 :: vp, vn ! flux-based cross-advective speeds (Yamazaki eq. 22)
+ real*4 :: qu
+ real*4 :: qd
+ real*4 :: un ! U_n advective speed (from east)
+ real*4 :: up ! U_p advective speed (from west)
+ real*4 :: umax ! local characteristic speed bound sqrt(g*hu) + |u|
+ real*4 :: dnminv ! 1/(D_nm + D_nmu) reused for both advective speeds
+ real*4 :: dzdx
+ !
+ real*4 :: hwet
+ real*4 :: phi
+ !
+ real*4 :: hu43
+ real*4 :: y_cbrt ! cube-root approximation hu^(1/3)
+ integer*4 :: i_cbrt ! bit pattern of hu/y_cbrt for the cube-root seed
+ !
+ real*4 :: zs2w, zs1e, dnm, dnmu, zrec ! advection work vars
+ real*4 :: zbup ! upwind still bed at u-point (bed of the cell the flow comes from)
+ integer :: ipw, ipe
+ real*4 :: zbnm, zbnmu ! bed at the west/east cell for subgrid advection (not necessarily zb)
+ real*4 :: mdrv ! subgrid wiggle suppression driver
+ !
+ real*4 :: min_dt_ip
+ !
+ logical :: iwet
+ !
+ call system_clock(count0, count_rate, count_max)
+ !
+ min_dt = dtmax
+ !
+ if (timestep_analysis) then
+ !
+ ! Do in loop for updating on GPU
+ !
+ !$acc parallel, present( timestep_analysis_required_timestep )
+ !$omp parallel &
+ !$omp private ( ip )
+ !$omp do
+ !$acc loop gang vector
+ do ip = 1, npuv
+ !
+ timestep_analysis_required_timestep(ip) = dtmax ! Reset per-cell limits; dry cells will retain dtmax
+ !
+ enddo
+ !$acc end parallel
+ !$omp end do
+ !$omp end parallel
+ !
+ endif
+ !
+ ! For some reason, it is necessary to set num_gangs here! Without, the program launches only 1 gang, and everything becomes VERY slow!
+ !
+ ! Copy velocity and flux from the previous time step (the velocity form
+ ! advects uv0; q0 provides the consistent previous-step fluxes for the
+ ! momentum-conserving cross-advection speeds)
+ !
+ !$acc parallel, present( uv, uv0, q, q0 )
+ !$omp parallel &
+ !$omp private ( ip )
+ !$omp do
+ !$acc loop gang vector
+ do ip = 1, npuv + ncuv
+ !
+ uv0(ip) = uv(ip)
+ q0(ip) = q(ip)
+ !
+ enddo
+ !$acc end parallel
+ !$omp end do
+ !$omp end parallel
+ !
+ !$omp parallel &
+ !$omp private ( ip,hu,ufr,nm,nmu,dzdx,frc,idir,itype,iref,dxuvinv,dxuv2inv,dyuvinv,dyuv2inv, &
+ !$omp uu_nm,uu_nmd,uu_nmu,uu_num,uu_ndm,vu, &
+ !$omp fcoriouv,gnavg2,iwet,zsu,dzuv,iuv,facint,fwmax,zmax,zmin,dqxudx,dqyudy,un,up,vp,vn,umax, &
+ !$omp dnminv,qu,qd,hwet,phi,adv,mdrv,hu43,y_cbrt,i_cbrt,min_dt_ip,zs2w,zs1e,dnm,dnmu,zrec,zbup,ipw,ipe,zbnm,zbnmu ) &
+ !$omp reduction ( min : min_dt )
+ !$omp do schedule ( dynamic, 256 )
+ !$acc parallel, present( kcuv, kfuv, zs, q, q0, uv, uv0, &
+ !$acc uv_flags_iref, uv_flags_type, uv_flags_dir, mask_adv, &
+ !$acc subgrid_uv_zmin, subgrid_uv_zmax, subgrid_uv_havg, subgrid_uv_nrep, subgrid_uv_pwet, &
+ !$acc subgrid_uv_havg_zmax, subgrid_uv_nrep_zmax, subgrid_uv_fnfit, subgrid_uv_navg_w, &
+ !$acc uv_index_z_nm, uv_index_z_nmu, uv_index_u_nmd, uv_index_u_nmu, uv_index_u_ndm, uv_index_u_num, &
+ !$acc uv_index_v_ndm, uv_index_v_ndmu, uv_index_v_nm, uv_index_v_nmu, cuv_index_uv, cuv_index_uv1, cuv_index_uv2, &
+ !$acc zb, tauwu, tauwv, patm, fwuv, gn2uv, dxminv, dxrinv, dyrinv, dxm2inv, dxr2inv, dyr2inv, &
+ !$acc dxrinvc, dyrinvc, fcorio2d, nuvisc, z_volume, cell_area, cell_area_m2, z_flags_iref, gnapp2, timestep_analysis_required_timestep ) num_gangs( 1024 ) vector_length( 128 )
+ !$acc loop, reduction( min : min_dt ), gang, vector
+ do ip = 1, npuv
+ !
+ if (kcuv(ip) == 1 .or. kcuv(ip) == 6) then
+ !
+ ! Regular UV point (or a coastal lateral boundary point)
+ !
+ ! Indices of surrounding water level points
+ !
+ nm = uv_index_z_nm(ip)
+ nmu = uv_index_z_nmu(ip)
+ !
+ iwet = .false.
+ !
+ ! Upwind surface at the u-point: take the surface from the cell the flow comes from
+ ! (sign of the previous-step velocity). The upwind bed (zbup) is only needed for the
+ ! regular-grid conveyance and is set in the non-subgrid branch below.
+ !
+ if (uv0(ip) > 1.0e-6) then
+ zsu = zs(nm)
+ elseif (uv0(ip) < -1.0e-6) then
+ zsu = zs(nmu)
+ else
+ zsu = max(zs(nm), zs(nmu))
+ endif
+ !
+ if (subgrid) then
+ !
+ zmin = subgrid_uv_zmin(ip)
+ zmax = subgrid_uv_zmax(ip)
+ !
+ if (zsu > zmin) then ! In the subgrid formulations, zmin is lowest pixel + huthresh. Huthresh was already applied when building the subgrid tables. In sfincs_domain, huthresh is set to 0.0 to acount for this.
+ iwet = .true.
+ endif
+ !
+ else
+ !
+ ! Flow depth at the u-point: D = upwind surface - upwind bed = zsu - zbup, i.e. the water
+ ! depth in the cell the flow comes from. The upwind bed follows the same flow direction
+ ! that selected zsu. The face is wet when that upwind depth exceeds huthresh. This avoids
+ ! the average-bed depth overshoot on steep downslopes while still allowing run-up.
+ !
+ if (uv0(ip) > 1.0e-6) then
+ zbup = zb(nm)
+ elseif (uv0(ip) < -1.0e-6) then
+ zbup = zb(nmu)
+ else
+ if (zs(nm) >= zs(nmu)) then
+ zbup = zb(nm)
+ else
+ zbup = zb(nmu)
+ endif
+ endif
+ !
+ if (zsu - zbup > huthresh) then
+ iwet = .true.
+ endif
+ !
+ endif
+ !
+ if (iwet) then
+ !
+ ! UV point is wet
+ !
+ if (use_quadtree) then
+ iref = uv_flags_iref(ip) ! refinement level
+ itype = uv_flags_type(ip) ! -1 is fine to coarse, 0 is normal, 1 is coarse to fine
+ else
+ iref = 1
+ itype = 0
+ endif
+ !
+ idir = uv_flags_dir(ip) ! 0 is u, 1 is v
+ !
+ ! Determine grid spacing (and coriolis factor fcoriouv)
+ !
+ if (crsgeo) then
+ !
+ ! Geographic coordinate system
+ !
+ if (itype==0) then
+ !
+ ! Regular
+ !
+ if (idir==0) then
+ !
+ ! U point
+ !
+ dxuvinv = dxminv(ip)
+ dyuvinv = dyrinv(iref)
+ dxuv2inv = dxm2inv(ip)
+ dyuv2inv = dyr2inv(iref)
+ !
+ else
+ !
+ ! V point
+ !
+ dxuvinv = dyrinv(iref)
+ dyuvinv = dxminv(ip)
+ dxuv2inv = dyr2inv(iref)
+ dyuv2inv = dxm2inv(ip)
+ !
+ endif
+ !
+ else
+ !
+ ! Fine to coarse or coarse to fine
+ !
+ if (idir==0) then
+ !
+ dxuvinv = 1.0 / (3*(1.0/dxminv(ip))/2)
+ dyuvinv = dyrinv(iref)
+ dxuv2inv = 0.0 ! no viscosity term
+ dyuv2inv = dyr2inv(iref)
+ !
+ else
+ !
+ dxuvinv = 1.0 / (3*(1.0/dyrinv(iref))/2)
+ dyuvinv = dxminv(ip)
+ dxuv2inv = 0.0 ! no viscosity term
+ dyuv2inv = dxm2inv(ip)
+ !
+ endif
+ !
+ endif
+ !
+ fcoriouv = fcorio2d(nm)
+ !
+ else
+ !
+ ! Projected coordinate system
+ !
+ if (itype==0) then
+ !
+ ! Regular
+ !
+ if (idir==0) then
+ !
+ ! U point
+ !
+ dxuvinv = dxrinv(iref)
+ dyuvinv = dyrinv(iref)
+ dxuv2inv = dxr2inv(iref)
+ dyuv2inv = dyr2inv(iref)
+ !
+ else
+ !
+ ! V point
+ !
+ dxuvinv = dyrinv(iref)
+ dyuvinv = dxrinv(iref)
+ dxuv2inv = dyr2inv(iref)
+ dyuv2inv = dxr2inv(iref)
+ !
+ endif
+ !
+ else
+ !
+ ! Fine to coarse or coarse to fine
+ !
+ if (idir==0) then
+ !
+ ! U point
+ !
+ dxuvinv = dxrinvc(iref)
+ dyuvinv = dyrinv(iref)
+ dxuv2inv = 0.0 ! no viscosity
+ dyuv2inv = dyr2inv(iref)
+ !
+ else
+ !
+ ! V point
+ !
+ dxuvinv = dyrinvc(iref)
+ dyuvinv = dxrinv(iref)
+ dxuv2inv = 0.0 ! no viscosity
+ dyuv2inv = dxr2inv(iref)
+ !
+ endif
+ !
+ endif
+ !
+ fcoriouv = fcorio
+ !
+ endif
+ !
+ ! Get velocities from the previous time step
+ !
+ if (advection .or. coriolis .or. viscosity .or. friction2d) then
+ !
+ ! Get the neighbors
+ !
+ uu_nm = uv0(ip)
+ uu_nmd = uv0(uv_index_u_nmd(ip))
+ uu_nmu = uv0(uv_index_u_nmu(ip))
+ uu_ndm = uv0(uv_index_u_ndm(ip))
+ uu_num = uv0(uv_index_u_num(ip))
+ vu = (uv0(uv_index_v_ndm(ip)) + uv0(uv_index_v_ndmu(ip)) + uv0(uv_index_v_nm(ip)) + uv0(uv_index_v_nmu(ip))) / 4
+ !
+ endif
+ !
+ ! Wet fraction phi (for non-subgrid or original subgrid approach phi should be 1.0)
+ !
+ phi = 1.0
+ !
+ ! Compute water depth at uv point
+ !
+ if (subgrid) then
+ !
+ if (zsu > zmax) then
+ !
+ ! Entire cell is wet, no interpolation from table needed for depth hu
+ !
+ hu = subgrid_uv_havg_zmax(ip) + zsu
+ !
+ if (wave_enhanced_roughness) then
+ !
+ ! Apparent roughness is computed in sfincs_wave_enhanced_roughness.f90. It is called by sfincs_bmi.f90.
+ ! Note: wave enhanced roughness is only done for uv points that are completely wet!
+ !
+ gnavg2 = gnapp2(ip)
+ !
+ else
+ !
+ ! Use fitting function for gnavg2
+ !
+ gnavg2 = subgrid_uv_navg_w(ip) - (subgrid_uv_navg_w(ip) - subgrid_uv_nrep_zmax(ip)) / (subgrid_uv_fnfit(ip) * (zsu - zmax) + 1.0)
+ !
+ endif
+ !
+ else
+ !
+ ! Interpolation required
+ !
+ dzuv = (zmax - zmin) / (subgrid_nlevels - 1) ! level size (is storing this in memory faster?)
+ iuv = min(int((zsu - zmin) / dzuv) + 1, subgrid_nlevels - 1) ! index of level below zsu
+ facint = (zsu - (zmin + (iuv - 1) * dzuv) ) / dzuv ! 1d interpolation coefficient
+ !
+ hu = subgrid_uv_havg(iuv, ip) + (subgrid_uv_havg(iuv + 1, ip) - subgrid_uv_havg(iuv, ip)) * facint ! grid-average depth
+ gnavg2 = subgrid_uv_nrep(iuv, ip) + (subgrid_uv_nrep(iuv + 1, ip) - subgrid_uv_nrep(iuv, ip)) * facint ! representative g*n^2
+ phi = subgrid_uv_pwet(iuv, ip) + (subgrid_uv_pwet(iuv + 1, ip) - subgrid_uv_pwet(iuv, ip)) * facint ! wet fraction
+ !
+ endif
+ !
+ else
+ !
+ hu = zsu - zbup ! Flow depth D = upwind zeta - upwind bed
+ gnavg2 = gn2uv(ip)
+ !
+ endif
+ !
+ ! Compute wet average depth hwet (used in wind and wave forcing)
+ !
+ hwet = hu / phi
+ !
+ ! FORCING TERMS
+ !
+ ! Pressure term
+ !
+ ! Apply slope limiter to dzdx (turned off by default)
+ !
+ if (slopelim < 9999.0) then
+ !
+ dzdx = min(max((zs(nmu) - zs(nm)) * dxuvinv, -slopelim), slopelim)
+ !
+ else
+ !
+ dzdx = (zs(nmu) - zs(nm)) * dxuvinv
+ !
+ endif
+ !
+ ! Velocity form: build frc directly as an acceleration [m/s^2]. Forces that scale
+ ! with depth (pressure, viscosity, Coriolis, atm) are written WITHOUT hu -- their hu
+ ! would only be divided out again. Only the surface stresses (wind, waves) keep a /hu.
+ !
+ frc = - g * dzdx
+ !
+ if (advection) then
+ !
+! if (mask_adv(ip) == 1) then
+ !
+ ! Momentum-conserved advection (Yamazaki, Kowalik & Cheung 2009,
+ ! eqs 18/20/22), VELOCITY form. Streamwise advective speeds from the Mader
+ ! upwind-zeta flux FLU = mean(U)*(upwind surface zeta + still-depth h), h=-zb;
+ ! zeta reconstructed 2nd-order (upwind face, +/-2 stencil) when co-directional,
+ ! else 1st-order. Cross term: flux-based two-sided upwind (eq. 22, see below).
+ ! 'adv' is a VELOCITY tendency [m/s^2] (added straight into frc).
+ !
+ ! One path for regular and subgrid bathymetry: on subgrid models zb
+ ! is the EFFECTIVE bed zs - z_volume/area maintained every step in
+ ! continuity (zb_effective; the file zb is not a valid conveyance
+ ! bed), so D = zs - zb is the subgrid cell-mean depth there. Then
+ ! zrec - zbnm = D_cell + 0.5*(upwind dzs) -- centered depth when
+ ! flat, with the upwind-surface boost at a front (subgrid-consistent
+ ! Mader reconstruction).
+ !
+ dnm = max(zs(nm) - zb(nm), 0.0)
+ dnmu = max(zs(nmu) - zb(nmu), 0.0)
+ !
+ zbnm = zb(nm)
+ zbnmu = zb(nmu)
+ !
+ ipw = uv_index_u_nmd(ip)
+ ipe = uv_index_u_nmu(ip)
+ !
+ ! Only use the +/-2 stencil when the neighbor is a REGULAR uv point
+ ! (<= npuv). At a quadtree refinement boundary uv_index_u_nm* points to a
+ ! COMBINED uv point (index > npuv), for which uv_index_z_* is out of bounds
+ ! (those arrays are sized npuv); fall back to 1st order there. Note ipw/ipe
+ ! are never 0 (sfincs_domain sets a missing neighbor to ip), so the old
+ ! "> 0" test never triggered the fallback.
+ !
+ if (ipw > 0 .and. ipw <= npuv) then
+ zs2w = zs(uv_index_z_nm(ipw))
+ else
+ zs2w = zs(nm)
+ endif
+ !
+ if (ipe > 0 .and. ipe <= npuv) then
+ zs1e = zs(uv_index_z_nmu(ipe))
+ else
+ zs1e = zs(nmu)
+ endif
+ !
+ if (uu_nmd >= 0.0) then
+ zrec = 0.5 * (zs2w + zs(nm))
+ else
+ zrec = zs(nm)
+ endif
+ !
+ qd = 0.5 * (uu_nmd + uu_nm) * max(zrec - zbnm, 0.0) ! FLU_p (west cell)
+ !
+ if (uu_nmu <= 0.0) then
+ zrec = 0.5 * (zs(nmu) + zs1e)
+ else
+ zrec = zs(nmu)
+ endif
+ !
+ qu = 0.5 * (uu_nm + uu_nmu) * max(zrec - zbnmu, 0.0) ! FLU_n (east cell)
+ !
+ dnm = max(dnm + dnmu, huthresh) ! D_nm + D_nmu
+ dnminv = 2.0 / dnm
+ !
+ ! Clamp the advective speeds to the local characteristic speed
+ ! sqrt(g*hu) + |u|, the bound the time-step limiter guarantees to
+ ! resolve. On subgrid grids the cell-mean depth D used in dnminv can
+ ! be far smaller than the conveyance depth hu, so q/D can otherwise
+ ! exceed the advective CFL and blow up (e.g. steep valleys).
+ !
+ umax = sqrt(g * hu) + abs(uv0(ip))
+ !
+ up = min(max(qd * dnminv, 0.0), umax) ! U_p (advective speed, from west)
+ un = max(min(qu * dnminv, 0.0), -umax) ! U_n (advective speed, from east)
+ !
+ dqxudx = ( up * (uu_nm - uu_nmd) + un * (uu_nmu - uu_nm) ) * dxuvinv
+ !
+ ! Cross-advection of u by v -- momentum-conserving two-sided upwind
+ ! (Yamazaki et al. 2009, eq. 22): advective speeds from the y-FLUXES
+ ! through the faces below/above the u-point (mean of the two flanking
+ ! v-point fluxes, previous time step q0), normalized by the same total
+ ! depth as the streamwise term (dnminv). Transport INTO the point from
+ ! either side contributes; at a collision line (v converging from both
+ ! sides) both terms stay active, where the velocity-average form below
+ ! gave ~zero cross-advection.
+ !
+ vp = min(max( 0.5 * (q0(uv_index_v_ndm(ip)) + q0(uv_index_v_ndmu(ip))) * dnminv, 0.0 ), umax)
+ vn = max(min( 0.5 * (q0(uv_index_v_nm(ip)) + q0(uv_index_v_nmu(ip))) * dnminv, 0.0 ), -umax)
+ !
+ dqyudy = ( vp * (uu_nm - uu_ndm) + vn * (uu_num - uu_nm) ) * dyuvinv
+ !
+ adv = - phi * (dqxudx + dqyudy) ! velocity tendency [m/s^2]
+ !
+! frc = frc + min(max(adv, -advlim), advlim) ! add limited advective acceleration
+ frc = frc + adv ! advective speeds are already clamped to umax above
+ !
+! endif
+ !
+ endif
+ !
+ ! Viscosity term
+ !
+ if (viscosity) then
+ !
+ if (itype == 0) then
+ !
+ frc = frc + nuvisc(iref) * ( (uu_nmu - 2*uu_nm + uu_nmd ) * dxuv2inv + (uu_num - 2*uu_nm + uu_ndm ) * dyuv2inv )
+ !
+ else
+ !
+ ! Increase viscosity to prevent instabilities on refinement (related to advection?)
+ !
+ frc = frc + nuviscfac * nuvisc(iref) * ( (uu_nmu - 2*uu_nm + uu_nmd ) * dxuv2inv + (uu_num - 2*uu_nm + uu_ndm ) * dyuv2inv )
+ !
+ endif
+ !
+ endif
+ !
+ ! Coriolis term
+ !
+ if (coriolis) then
+ !
+ if (idir==0) then
+ !
+ frc = frc + fcoriouv * vu ! U
+ !
+ else
+ !
+ frc = frc - fcoriouv * vu ! V
+ !
+ endif
+ !
+ endif
+ !
+ ! Wind forcing
+ !
+ if (wind) then
+ !
+ if (hwet > 0.25) then
+ !
+ ! Wind stress is a surface force -> acceleration = stress / depth
+ !
+ if (idir==0) then
+ !
+ frc = frc + phi * tauwu(nm) / max(hu, huvmin)
+ !
+ else
+ !
+ frc = frc + phi * tauwv(nm) / max(hu, huvmin)
+ !
+ endif
+ !
+ else
+ !
+ ! Reduce wind drag at water depths < 0.25 m (tauw*hu*4 / hu = tauw*4)
+ !
+ if (idir==0) then
+ !
+ frc = frc + tauwu(nm) * 4
+ !
+ else
+ !
+ frc = frc + tauwv(nm) * 4
+ !
+ endif
+ !
+ endif
+ !
+ endif
+ !
+ ! Atmospheric pressure
+ !
+ if (patmos) then
+ !
+ frc = frc + (patm(nm) - patm(nmu)) * dxuvinv / rhow
+ !
+ endif
+ !
+ ! Wave forcing
+ !
+ if (snapwave) then
+ !
+ ! Limited wave forces in shallow water
+ !
+ ! facmax = 0.25*sqrt(g)*rhow*gammax**2
+ ! fmax = facmax*hu*sqrt(hu)/tp/rhow (we already divided by rhow in sfincs_snapwave)
+ !
+ fwmax = 0.8 * hwet * sqrt(hwet) / 15
+ !
+ ! Wave force is a surface force -> acceleration = force / depth
+ !
+ frc = frc + phi * sign(min(abs(fwuv(ip)), fwmax), fwuv(ip)) / max(hu, huvmin)
+ !
+ endif
+ !
+ ! hu**(1/3) and hu**(4/3) for the velocity-form Manning friction, via a fast cube
+ ! root: an integer bit-hack seed (magic constant 709921077) refined by one Newton
+ ! iteration (y <- y - (y^3 - hu)/(3 y^2)). ~0.1% accurate over the depth range,
+ ! no pow, no table. y_cbrt = hu^(1/3) is reused for the newly-wet estimate below.
+ !
+ i_cbrt = transfer(hu, i_cbrt)
+ i_cbrt = i_cbrt / 3 + 709921077
+ y_cbrt = transfer(i_cbrt, y_cbrt)
+ y_cbrt = y_cbrt - (y_cbrt * y_cbrt * y_cbrt - hu) / (3.0 * y_cbrt * y_cbrt)
+ hu43 = hu * y_cbrt
+ !
+ ! Friction velocity proxy ufr (velocity form: the implicit Manning factor is
+ ! gnavg2*ufr/hu^(4/3) with ufr the friction-driving velocity magnitude).
+ !
+ if (kfuv(ip) == 0) then
+ !
+ ! This uv point just became wet, so estimate the equilibrium velocity
+ ! (hu^(2/3) = (hu^(1/3))^2 = y_cbrt^2, reusing the cube root above)
+ !
+ ufr = sqrt(abs(dzdx) / (max(gnavg2, 1.0e-5) / 10)) * y_cbrt * y_cbrt
+ !
+ else
+ !
+ if (friction2d) then
+ !
+ ! Both velocity components: ufr = sqrt(u^2 + v^2)
+ !
+ ufr = sqrt(uv0(ip)**2 + vu**2)
+ !
+ else
+ !
+ ! Streamwise velocity only: ufr = |u|
+ !
+ ufr = abs(uv0(ip))
+ !
+ endif
+ !
+ endif
+ !
+ ! Velocity update. frc is a velocity tendency and the
+ ! implicit friction factor is the velocity-form Manning term gnavg2*|u|/hu^(4/3).
+ !
+ uv(ip) = (uv0(ip) + frc * dt) / (1.0 + gnavg2 * dt * ufr / hu43)
+ !
+ if (subgrid .and. wiggle_suppression) then
+ !
+ ! If the acceleration of water level in cell nm is large and positive and in nmu large and negative, or vice versa, apply limiter to the flux. Only for subgrid.
+ !
+ mdrv = abs(zsderv(nm) - zsderv(nmu)) - wiggle_threshold
+ !
+ if (mdrv > 0.0) then
+ !
+ uv(ip) = uv(ip) * wiggle_threshold / (wiggle_factor * mdrv + wiggle_threshold)
+ !
+ endif
+ !
+ endif
+ !
+ ! Velocity limiter (default 10 m/s)
+ !
+ uv(ip) = min(max(uv(ip), - uvlim), uvlim)
+ !
+ ! No flow out of a cell that is (going) dry
+ !
+ if (zs(nm) < zb(nm)) uv(ip) = min(uv(ip), 0.0)
+ if (zs(nmu) < zb(nmu)) uv(ip) = max(uv(ip), 0.0)
+ !
+ ! Continuity flux from the updated velocity and the conveyance depth.
+ !
+ q(ip) = uv(ip) * hu
+ !
+ kfuv(ip) = 1
+ !
+ ! Determine minimum time step (alpha is added later on in sfincs_lib.f90) of all uv points
+ ! Use maximum of sqrt(gh) and current velocity
+ !
+ min_dt_ip = 1.0 / ( max(sqrt(g * hu), abs(uv(ip)) ) * dxuvinv)
+ !
+ min_dt = min(min_dt, min_dt_ip)
+ !
+ ! Compute timestep per grid cell
+ !
+ if (timestep_analysis) then
+ !
+ timestep_analysis_required_timestep(ip) = min_dt_ip
+ !
+ endif
+ !
+ else
+ !
+ q(ip) = 0.0
+ uv(ip) = 0.0
+ kfuv(ip) = 0
+ !
+ endif
+ !
+ endif
+ enddo
+ !$omp end do
+ !$omp end parallel
+ !$acc end parallel
+ !
+ if (ncuv > 0) then
+ !
+ ! Loop through combined uv points and determine average uv and q
+ ! The combined q and uv values are used in the continuity equation and in the netcdf output
+ !
+ !$omp parallel &
+ !$omp private ( icuv )
+ !$omp do
+ !$acc parallel, present( q, uv, cuv_index_uv, cuv_index_uv1, cuv_index_uv2 )
+ !$acc loop gang vector
+ do icuv = 1, ncuv
+ !
+ ! Average of the two uv points
+ !
+ q(cuv_index_uv(icuv)) = (q(cuv_index_uv1(icuv)) + q(cuv_index_uv2(icuv))) / 2
+ uv(cuv_index_uv(icuv)) = (uv(cuv_index_uv1(icuv)) + uv(cuv_index_uv2(icuv))) / 2
+ !
+ enddo
+ !$acc end parallel
+ !$omp end do
+ !$omp end parallel
+ !
+ endif
+ !
+ call system_clock(count1, count_rate, count_max)
+ tloop = tloop + 1.0*(count1 - count0)/count_rate
+ !
+ end subroutine
+ !
+end module
diff --git a/source/src/sfincs_ncoutput.F90 b/source/src/sfincs_ncoutput.F90
index 04427c938..fd3dc9c39 100644
--- a/source/src/sfincs_ncoutput.F90
+++ b/source/src/sfincs_ncoutput.F90
@@ -1702,11 +1702,11 @@ subroutine ncoutput_add_params(ncid, varid)
NF90(nf90_put_att(ncid, varid, 'factor_pres',factor_pres))
NF90(nf90_put_att(ncid, varid, 'factor_prcp',factor_prcp))
NF90(nf90_put_att(ncid, varid, 'factor_spw_size',factor_spw_size))
- NF90(nf90_put_att(ncid, varid, 'nonh',logical2int(nonhydrostatic)))
- NF90(nf90_put_att(ncid, varid, 'nh_fnudge',nh_fnudge))
- NF90(nf90_put_att(ncid, varid, 'nh_tstop',nh_tstop))
- NF90(nf90_put_att(ncid, varid, 'nh_tol',nh_tol))
- NF90(nf90_put_att(ncid, varid, 'nh_itermax',nh_itermax))
+ !NF90(nf90_put_att(ncid, varid, 'nonh',logical2int(nonhydrostatic)))
+ !NF90(nf90_put_att(ncid, varid, 'nonh_fnudge',nonh_fnudge))
+ !NF90(nf90_put_att(ncid, varid, 'nonh_tstop',nonh_tstop))
+ !NF90(nf90_put_att(ncid, varid, 'nonh_tol',nonh_tol))
+ !NF90(nf90_put_att(ncid, varid, 'nonh_itermax',nonh_itermax))
!
! Domain
!
diff --git a/source/src/sfincs_nonhydrostatic.f90 b/source/src/sfincs_nonhydrostatic.f90
index 6200a1b5b..8237a9ee0 100644
--- a/source/src/sfincs_nonhydrostatic.f90
+++ b/source/src/sfincs_nonhydrostatic.f90
@@ -1,49 +1,189 @@
-! Non-hydrostatic code now only works with regular grids (can still use quadtree netcdf file as long as there are no refinement levels).
-! Now uses bicgstab_ilu to solve matrix. Both should ideally utilize CPU and GPU parallelization. Currently, this solver cannot be fully parallelized.
-! Should try to find a solver that can be.
-
+! Non-hydrostatic pressure projection (single-layer Keller-box), GPU-ready.
+!
+! The depth-averaged momentum/continuity set is closed with a linearly-varying
+! non-hydrostatic pressure pnh (bed value; depth mean pnh/2). Eliminating the
+! provisional velocities gives a symmetric-positive-definite pressure operator
+!
+! A = (dt/rho) G^T diag(hu) G + diag( kbfac dt / (rho H) ) ,
+!
+! built on the momentum stencil (G and the divergence are exact negative
+! transposes), solved with a Chronopoulos-Gear (single-reduction) conjugate
+! gradient, warm-started from the previous step. Works on regular and quadtree
+! grids (refinement-transition cells are excluded from the nonh mask and stay
+! hydrostatic, so every nonh face joins two same-level cells) and on projected
+! and geographic (crsgeo) grids via per-face/per-row metric spacing.
+!
+! The solver machinery is designed for memory bandwidth and GPU offload:
+!
+! 1) ASSEMBLED 5-POINT STENCIL. The symmetric operator is assembled ONCE PER
+! SOLVE into four off-diagonal arrays (nh_aw/ae/as/an) with STATIC neighbour
+! indices (nh_nbw/nbe/nbs/nbn). The matvec is then a single loop:
+! w(i) = r(i) + aw*r(nbw) + ae*r(nbe) + as*r(nbs) + an*r(nbn)
+! instead of a per-iteration two-pass face form.
+!
+! 2) SYMMETRIC JACOBI SCALING BAKED IN. The system is solved as
+! (S A S) y = S b, pnh = S y, S = diag(1/sqrt(diag(A))),
+! so the scaled operator has UNIT diagonal: no preconditioner work in the
+! iteration, and = removes one reduction term.
+!
+! 3) BREAKING DIRICHLET BAKED IN. The row AND column of a breaking cell are
+! zeroed at assembly time (nh_sclo = nh_scl * activity, symmetric), the
+! scaled diagonal stays 1, and b=0 there -> pnh=0 is held exactly with NO
+! per-iteration branch (faces stay active for the flux correction).
+!
+! 4) GHOST ELEMENT 0. CG vectors run (0:nrows) with element 0 pinned at zero;
+! a missing neighbour has index 0 and a zero coefficient (nh_sclo(0)=0), so
+! the matvec has NO conditionals at all -- fully vectorizable on CPU,
+! coalesced on GPU.
+!
+! 5) TWO FUSED KERNELS PER ITERATION:
+! kernel B: w = A r fused with the dots gam=, del=
+! kernel A: p=r+beta*p ; s=w+beta*s ; y+=alpha*p ; r-=alpha*s
+! One host synchronisation point per iteration (the reduction). Dot products
+! accumulate in real*8 (the real*4 vectors are unchanged).
+!
+! 6) DUAL DIRECTIVES. Every kernel carries both !$omp (CPU, schedule(static))
+! and !$acc (GPU) directives, following the sfincs_momentum.f90 convention.
+! All module arrays are made device-resident at initialization (enter data);
+! per-step work then runs entirely on the device. The host copy of pnh is
+! refreshed only at map output times ('update host(pnh) if_present' in
+! sfincs_output, alongside the other output transfers).
+!
+! NOTE: the convergence test uses the SCALED residual norm ||S r|| / ||S b||
+! (the natural Jacobi-CG norm); nonh_tol applies to that.
+!
+! On subgrid models zb is the EFFECTIVE bed zs - z_volume/area maintained every
+! step in continuity (zb_effective), so all depth checks here see the volume-
+! consistent depth. The bed SLOPES for the bottom kinematic condition w_b are
+! FROZEN at initialization (nh_slbed): file bed on regular grids, subgrid
+! cell-mean bed z_zmax - z_volmax/area on subgrid models.
+!
module sfincs_nonhydrostatic
- !
- integer, dimension(:,:), allocatable :: index_sparse_matrix
- integer, dimension(:,:), allocatable :: nh_uv_index
+ !
+ implicit none
+ !
+ ! Input parameters (sfincs.inp, read in sfincs_input)
+ !
+ real*4 :: nonh_fnudge ! fraction of the pressure correction applied to uv/q
+ real*4 :: nonh_tstop ! time (s, absolute after input processing) to stop applying nonh corrections
+ integer :: nonh_itermax ! max CG iterations per solve
+ integer :: nonh_fadein ! open-boundary nonh fade-in width (cells); 0 = off
+ real*4 :: nonh_filter ! global 2dx pnh filter weight; 0 = off
+ real*4 :: nonh_dzbmax ! cap on |d(zb)/dx| in the bottom kinematic w_b; 0 = no cap
+ real*4 :: nonh_brsteep ! HFA steepness breaking onset; 0 = off
+ real*4 :: nonh_brfr ! Froude breaking onset (NEOWAVE); 0 = off -> use nonh_brsteep
+ integer :: nonh_brsmooth ! breaking-flag smoothing passes: ramps pnh out over ~brsmooth+1 cells at the breaking-zone edges; 0 = sharp (one face)
+ integer :: nonh_slsmooth ! frozen-bed-slope smoothing passes: bounds bed curvature d2zb/dx2 at slope breaks (e.g. island toe); 0 = off
+ real*4 :: nonh_treform ! breaking reformation time scale (s): released cells recover gradually, brfac += dt/treform; 0 = instant
+ real*4 :: nonh_smoothbnd ! localized 2dx pnh smoothing strength (fade-in / shallow zones)
+ real*4 :: nonh_smoothdep ! depth below which the localized smoothing also acts; 0 = off
+ real*4 :: nonh_disp ! Keller-box vertical factor (dispersion tuning; 1.0 ~ Airy, 2.0 strict)
+ real*4 :: nonh_pmax ! depth limiter: |pnh| <= nonh_pmax * rho g H; 0 = off
+ real*4 :: nonh_tol ! CG relative tolerance (scaled residual norm)
+ logical :: nonh_movingbed ! add d(zb)/dt (= dzbext/dt) to the bottom kinematic w_b, so a moving seafloor radiates a depth-filtered (Kajiura-like) surface response
+ !
+ ! Non-hydrostatic cell mask (filled from the quadtree file in sfincs_domain)
+ !
+ integer*1, dimension(:), allocatable :: mask_nonh
+ !
integer, dimension(:), allocatable :: nm_index_of_row
integer, dimension(:), allocatable :: row_index_of_nm
integer, dimension(:), allocatable :: uv_index_of_nhuv
- integer, dimension(:), allocatable :: col_idx
- integer, dimension(:), allocatable :: row_ptr
!
- real*4, dimension(:), allocatable :: pnh
- real*4, dimension(:), allocatable :: ws
- real*4, dimension(:), allocatable :: wb
- real*4, dimension(:), allocatable :: wb0
+ real*4, dimension(:), allocatable :: pnh ! non-hydrostatic bed pressure (nrows)
+ real*4, dimension(:), allocatable :: ws ! surface vertical velocity (nrows)
+ real*4, dimension(:), allocatable :: wb ! bottom vertical velocity (nrows)
+ real*4, dimension(:), allocatable :: wb0 ! previous bottom v. velocity (nrows)
+ !
+ real*4, dimension(:), allocatable :: Dnm ! layer depth per row (nrows)
+ !
+ ! Per-row inverse grid spacing in metres (handles quadtree refinement levels and
+ ! geographic crsgeo lat-dependence). nh_dxr/nh_dyr are the per-CELL 1/dx, 1/dy used
+ ! by the breaking dzdt divergence and the frozen bed slopes; nh_cf (below) is the
+ ! per-FACE 0.5/dx used in the pressure operator.
+ !
+ real*4, dimension(:), allocatable :: nh_dxr ! per-row 1/dx in metres (nrows)
+ real*4, dimension(:), allocatable :: nh_dyr ! per-row 1/dy in metres (nrows)
+ !
+ ! Bed slopes for the bottom kinematic condition w_b (step 7), FROZEN at
+ ! initialization (file bed on regular grids, subgrid cell-mean bed
+ ! z_zmax - z_volmax/area on subgrid models) and capped at +/- nonh_dzbmax.
+ ! They must NOT follow the per-step effective bed zs - z_volume/area that
+ ! subgrid models maintain (zb_effective): that bed moves with the wet fraction
+ ! every step and its slope jitter would feed straight into w_b/ws. The real
+ ! bed slope is constant in time.
+ !
+ real*4, dimension(:,:), allocatable :: nh_slbed ! (4, nrows) slots: 1 left(md), 2 right(mu), 3 below(nd), 4 above(nu)
+ !
+ real*4, dimension(:), allocatable :: nh_dvert ! vertical-accel diagonal term (nrows)
+ real*4, dimension(:), allocatable :: nh_fade ! open-boundary nonh fade-in 0..1 (nrows)
+ real*4, dimension(:), allocatable :: nh_brfac ! breaking STATE 1=full nonh .. 0=hydrostatic (nrows; hysteresis memory, unsmoothed)
+ real*4, dimension(:), allocatable :: nh_bract ! breaking ACTIVITY used by the pressure operator (nrows; = nh_brfac, optionally smoothed)
+ !
+ integer, dimension(:), allocatable :: nh_faceuv ! full uv index of each nh face (nhuv)
+ integer, dimension(:), allocatable :: nh_faceL ! row index of left/bottom cell (nhuv, 0 = boundary)
+ integer, dimension(:), allocatable :: nh_faceR ! row index of right/top cell (nhuv, 0 = boundary)
+ real*4, dimension(:), allocatable :: nh_cf ! 0.5/dx static gradient weight (nhuv)
+ real*4, dimension(:), allocatable :: nh_cR ! gradient coeff of p(R) per face(nhuv)
+ real*4, dimension(:), allocatable :: nh_cL ! gradient coeff of p(L) per face(nhuv)
+ real*4, dimension(:), allocatable :: nh_hu ! water depth at face this step (nhuv)
+ !
+ ! Per-row incident-face map. For each row, the (up to 4) non-hydrostatic faces
+ ! touching it, and implicitly which side the row is on:
+ ! slot 1 = left face -> row is the R-cell -> use nh_cR ; "other" cell = nh_faceL
+ ! slot 2 = right face -> row is the L-cell -> use nh_cL ; "other" cell = nh_faceR
+ ! slot 3 = below face -> row is the R-cell -> use nh_cR ; "other" cell = nh_faceL
+ ! slot 4 = above face -> row is the L-cell -> use nh_cL ; "other" cell = nh_faceR
+ ! 0 = no face on that side.
+ !
+ integer, dimension(:,:), allocatable :: nh_cellface ! (4, nrows)
+ !
+ ! Static 5-point neighbour map: row index of the West/East/South/North
+ ! neighbour of each row (0 = no nonh neighbour on that side -> ghost row 0).
+ !
+ integer, dimension(:), allocatable :: nh_nbw, nh_nbe, nh_nbs, nh_nbn ! (nrows)
+ !
+ ! Assembled, symmetrically-scaled off-diagonals of the operator (unit diagonal),
+ ! rebuilt every solve (wet/dry, depths and breaking change each step).
!
- real*4, dimension(:), allocatable :: Dnm
- real*4, dimension(:), allocatable :: dzbdx
- real*4, dimension(:), allocatable :: dzbdy
+ real*4, dimension(:), allocatable :: nh_aw, nh_ae, nh_as, nh_an ! (nrows)
+ real*4, dimension(:), allocatable :: nh_scl ! 1/sqrt(diag(A)) (nrows)
+ real*4, dimension(:), allocatable :: nh_sclo ! scl * breaking-activity (0:nrows), ghost 0
+ !
+ ! CG vectors, ghost element 0 pinned at zero (never written by any kernel).
+ !
+ real*4, dimension(:), allocatable :: nh_y ! scaled pressure iterate (0:nrows)
+ real*4, dimension(:), allocatable :: nh_r ! residual (0:nrows)
+ real*4, dimension(:), allocatable :: nh_p ! search direction (0:nrows)
+ real*4, dimension(:), allocatable :: nh_s ! s = A p (recurrence) (0:nrows)
+ real*4, dimension(:), allocatable :: nh_w ! w = A r (0:nrows)
+ real*4, dimension(:), allocatable :: nh_b ! scaled right-hand side (0:nrows)
!
real*4 :: huthresh_nh
!
integer :: nrows
- integer :: nr_vals_in_matrix
integer :: nhuv
!
+ ! CG iteration diagnostics (cumulative over the run; reported in the timing summary)
+ !
+ integer*8 :: nh_iter_total = 0 ! sum of CG iterations over all solves
+ integer*8 :: nh_solve_count = 0 ! number of CG solves (active time steps)
+ integer :: nh_iter_max = 0 ! most CG iterations in any single solve
+ !
contains
!
subroutine initialize_nonhydrostatic()
!
- ! Initialization of pardiso solver
+ ! Builds the row<->cell maps, the nonh face maps, the static neighbour map,
+ ! the frozen bed slopes and the open-boundary fade-in, allocates the solver
+ ! work arrays and pushes everything the per-step kernels touch onto the GPU.
!
use sfincs_data
+ use sfincs_log
!
implicit none
!
- integer nm, k, nmu, nmd, num, ndm, irow, icol, inb, ip
- integer inhuv, iii
- !
- ! Temporary arrays
- !
- integer, dimension(:), allocatable :: col_idx0
- integer, dimension(:,:), allocatable :: nh_nm_index
+ integer :: nm, irow, ip, inhuv, rl, rr, i, f
!
allocate(row_index_of_nm(np))
!
@@ -67,37 +207,28 @@ subroutine initialize_nonhydrostatic()
allocate(ws(nrows))
allocate(wb(nrows))
allocate(wb0(nrows))
- allocate(index_sparse_matrix(5, nrows))
- allocate(row_ptr(nrows + 1))
- allocate(col_idx0(5 * nrows))
allocate(Dnm(nrows))
- allocate(dzbdx(nrows))
- allocate(dzbdy(nrows))
allocate(nm_index_of_row(nrows))
- allocate(nh_uv_index(4, nrows))
- allocate(nh_nm_index(4, nrows))
+ allocate(nh_dxr(nrows))
+ allocate(nh_dyr(nrows))
+ allocate(nh_dvert(nrows))
+ allocate(nh_fade(nrows))
+ allocate(nh_brfac(nrows))
+ allocate(nh_bract(nrows))
+ !
+ pnh = 0.0
+ ws = 0.0
+ wb = 0.0
+ wb0 = 0.0
+ Dnm = 0.0
+ nh_dvert = 0.0
+ nh_fade = 1.0
+ nh_brfac = 1.0
+ nh_bract = 1.0
+ nm_index_of_row = 0
!
huthresh_nh = max(huthresh, 0.01)
!
- pnh = 0.0
- ws = 0.0
- wb = 0.0
- wb0 = 0.0
- !
- index_sparse_matrix = 0
- row_ptr = 0
- col_idx = 0
- col_idx0 = 0
- irow = 0
- k = 0
- Dnm = 0.0
- dzbdx = 0.0
- dzbdy = 0.0
- nm_index_of_row = 0
- nh_uv_index = 0
- nh_nm_index = 0
- ! bnd = 0
- !
! Map row indices to nm indices and vice versa
!
irow = 0
@@ -107,24 +238,22 @@ subroutine initialize_nonhydrostatic()
irow = irow + 1
row_index_of_nm(nm) = irow
nm_index_of_row(irow) = nm
+ !
+ ! Per-cell inverse grid spacing in metres. y-spacing is uniform per refinement
+ ! level (dyrinv) in both projected and geographic mode. x-spacing is uniform per
+ ! level (dxrinv) when projected, but lat-dependent per cell (1/dxm) when crsgeo.
+ !
+ nh_dyr(irow) = dyrinv(z_flags_iref(nm))
+ if (crsgeo) then
+ nh_dxr(irow) = 1.0 / dxm(nm)
+ else
+ nh_dxr(irow) = dxrinv(z_flags_iref(nm))
+ endif
endif
- enddo
- !
- ! Find velocity points needed for nh computations
- !
- ! 4
- ! +------|------+
- ! | |
- ! | irow |
- ! 1 - + 5 - 2
- ! | |
- ! | |
- ! +------|------+
- ! 3
- !
- ! Loop through all uv points to get nh uv points (any uv point that touches nh cell)
+ enddo
!
- ! First just count them
+ ! Find velocity points needed for nh computations: any uv point that touches
+ ! a nonh cell. First count, then store.
!
nhuv = 0
!
@@ -140,607 +269,953 @@ subroutine initialize_nonhydrostatic()
!
allocate(uv_index_of_nhuv(nhuv))
uv_index_of_nhuv = 0
- !
- ! Now find neigboring indices of nh uv points
!
inhuv = 0
!
do ip = 1, npuv
!
- nm = uv_index_z_nm(ip)
- nmu = uv_index_z_nmu(ip)
- !
- if (mask_nonh(nm) == 1 .or. mask_nonh(nmu) == 1) then
+ if (mask_nonh(uv_index_z_nm(ip)) == 1 .or. mask_nonh(uv_index_z_nmu(ip)) == 1) then
!
inhuv = inhuv + 1
- !
- uv_index_of_nhuv(inhuv) = ip ! uv index of this nhuv point
- !
- ! Get direction of point
- !
- if (uv_flags_dir(ip) == 0) then
- !
- ! x
- !
- if (mask_nonh(nm) == 1) then
- !
- ! nm point
- !
- irow = row_index_of_nm(nm)
- inb = row_index_of_nm(nmu)
- !
- nh_uv_index(2, irow) = inhuv
- nh_nm_index(2, irow) = inb
- !
- endif
- !
- if (mask_nonh(nmu) == 1) then
- !
- ! nmu point
- !
- inb = row_index_of_nm(nm)
- irow = row_index_of_nm(nmu)
- !
- nh_uv_index(1, irow) = inhuv
- nh_nm_index(1, irow) = inb
- !
- endif
- !
- else
- !
- ! y
- !
- if (mask_nonh(nm) == 1) then
- !
- ! nm point
- !
- irow = row_index_of_nm(nm)
- inb = row_index_of_nm(nmu)
- !
- nh_uv_index(4, irow) = inhuv
- nh_nm_index(4, irow) = inb
- !
- endif
- !
- if (mask_nonh(nmu) == 1) then
- !
- ! nmu point
- !
- inb = row_index_of_nm(nm)
- irow = row_index_of_nm(nmu)
- !
- nh_uv_index(3, irow) = inhuv
- nh_nm_index(3, irow) = inb
- !
- endif
- !
- endif
+ uv_index_of_nhuv(inhuv) = ip
!
endif
!
- enddo
+ enddo
!
- ! Determine indices of sparse matrix AA
+ ! Face arrays
+ !
+ allocate(nh_faceuv(nhuv))
+ allocate(nh_faceL(nhuv))
+ allocate(nh_faceR(nhuv))
+ allocate(nh_cf(nhuv))
+ allocate(nh_cR(nhuv))
+ allocate(nh_cL(nhuv))
+ allocate(nh_hu(nhuv))
+ nh_faceuv = 0
+ nh_faceL = 0
+ nh_faceR = 0
+ nh_cf = 0.0
+ nh_cR = 0.0
+ nh_cL = 0.0
+ nh_hu = 0.0
+ !
+ ! Static per-face data. nh_cf = 0.5 / dx is the (depth-averaged) gradient
+ ! weight; the factor 0.5 reflects that the linearly-varying non-hydrostatic
+ ! pressure has a depth mean of p_bed/2. The spacing is taken PER FACE so it is
+ ! correct across quadtree refinement levels and for geographic (crsgeo) grids:
+ ! x-faces use the per-uv-point metre spacing dxminv (crsgeo) or the per-level
+ ! dxrinv (projected), y-faces the per-level dyrinv (uniform in both modes).
+ ! nh_cf is a single value per face, used identically by the gradient and its
+ ! transpose divergence, so the operator A = G^T diag(hu) G stays symmetric
+ ! regardless of spacing variation.
!
- iii = 0
- do irow = 1, nrows
- !
- ! if (irow <= iii) then
- ! bnd(icol) = 1
- ! endif
- !
- nm = nm_index_of_row(irow)
- !
- ! Left
- !
- icol = nh_nm_index(1, irow)
- !
- if (icol > iii) then
- !
- ! Cell has a neighbor to the left
- !
- k = k + 1
- col_idx0(k) = icol
- index_sparse_matrix(1, irow) = k
- !
- if (row_ptr(irow) == 0) then ! first data in row
- !
- row_ptr(irow) = k
- !
- endif
- !
- endif
- !
- ! Bottom
- !
- icol = nh_nm_index(3, irow)
- !
- if (icol > iii) then
- !
- ! Cell has a neighbor below
- !
- k = k + 1
- col_idx0(k) = icol
- index_sparse_matrix(3, irow) = k
- !
- if (row_ptr(irow) == 0) then ! first data in row
- !
- row_ptr(irow) = k
- !
- endif
- !
- endif
- !
- ! Centre
- !
- icol = irow
- !
- ! if (icol > iii) then
- !
- k = k + 1
- col_idx0(k) = icol
- index_sparse_matrix(5, irow) = k
- !
- if (row_ptr(irow) == 0) then ! first data in row
- !
- row_ptr(irow) = k
- !
- endif
- !
- ! endif
- !
- ! Top
+ do ip = 1, nhuv
!
- icol = nh_nm_index(4, irow)
+ inhuv = uv_index_of_nhuv(ip)
+ nh_faceuv(ip) = inhuv
+ nh_faceL(ip) = row_index_of_nm(uv_index_z_nm(inhuv))
+ nh_faceR(ip) = row_index_of_nm(uv_index_z_nmu(inhuv))
!
- if (icol > iii) then
- !
- ! Cell has a neighbor above
- !
- k = k + 1
- col_idx0(k) = icol
- index_sparse_matrix(4, irow) = k
- !
- if (row_ptr(irow) == 0) then ! first data in row
- !
- row_ptr(irow) = k
- !
- endif
- !
+ if (uv_flags_dir(inhuv) == 0) then
+ if (crsgeo) then
+ nh_cf(ip) = 0.5 * dxminv(inhuv)
+ else
+ nh_cf(ip) = 0.5 * dxrinv(uv_flags_iref(inhuv))
+ endif
+ else
+ nh_cf(ip) = 0.5 * dyrinv(uv_flags_iref(inhuv))
endif
!
- ! Right
- !
- icol = nh_nm_index(2, irow)
- !
- if (icol > iii) then
- !
- ! Cell has a neighbor to the right
- !
- k = k + 1
- col_idx0(k) = icol
- index_sparse_matrix(2, irow) = k
- !
- if (row_ptr(irow) == 0) then ! first data in row
- !
- row_ptr(irow) = k
- !
- endif
- !
+ enddo
+ !
+ ! Per-row incident-face map (see header). Walk the faces once: each face is
+ ! the RIGHT face of its L-cell and the LEFT face of its R-cell (x), or the
+ ! ABOVE face of its L-cell and the BELOW face of its R-cell (y).
+ !
+ allocate(nh_cellface(4, nrows))
+ nh_cellface = 0
+ do ip = 1, nhuv
+ inhuv = nh_faceuv(ip)
+ rl = nh_faceL(ip)
+ rr = nh_faceR(ip)
+ if (uv_flags_dir(inhuv) == 0) then
+ if (rl > 0) nh_cellface(2, rl) = ip ! face is the L-cell's RIGHT face
+ if (rr > 0) nh_cellface(1, rr) = ip ! face is the R-cell's LEFT face
+ else
+ if (rl > 0) nh_cellface(4, rl) = ip ! face is the L-cell's ABOVE face
+ if (rr > 0) nh_cellface(3, rr) = ip ! face is the R-cell's BELOW face
endif
- !
enddo
!
- nr_vals_in_matrix = k
- allocate(col_idx(nr_vals_in_matrix))
- col_idx(1:k) = col_idx0(1:k)
- row_ptr(nrows + 1) = k + 1
+ ! Static 5-point neighbour map: neighbour row of each row through its (up to 4)
+ ! incident nonh faces. nh_faceL/R are 0 on a domain/open-boundary side ->
+ ! neighbour 0 -> ghost element.
!
- ! Compute bed level slopes
+ allocate(nh_nbw(nrows))
+ allocate(nh_nbe(nrows))
+ allocate(nh_nbs(nrows))
+ allocate(nh_nbn(nrows))
!
- do irow = 1, nrows
- !
- nm = nm_index_of_row(irow)
- !
- nmd = 0
- nmu = 0
- ndm = 0
- num = 0
- !
- ! Left
- !
- if (nh_nm_index(1, irow) > 0) then
- !
- nmd = nm_index_of_row(nh_nm_index(1, irow))
- dzbdx(irow) = dzbdx(irow) + 0.5 * (zb(nm) - zb(nmd)) * dxrinv(1)
- !
- endif
+ do i = 1, nrows
+ nh_nbw(i) = 0 ; f = nh_cellface(1, i) ; if (f > 0) nh_nbw(i) = nh_faceL(f)
+ nh_nbe(i) = 0 ; f = nh_cellface(2, i) ; if (f > 0) nh_nbe(i) = nh_faceR(f)
+ nh_nbs(i) = 0 ; f = nh_cellface(3, i) ; if (f > 0) nh_nbs(i) = nh_faceL(f)
+ nh_nbn(i) = 0 ; f = nh_cellface(4, i) ; if (f > 0) nh_nbn(i) = nh_faceR(f)
+ enddo
+ !
+ ! Frozen, capped per-face bed slopes for the bottom kinematic condition
+ ! w_b = u . d(zb)/dx (step 7). Computed here from the initialization-time zb
+ ! (file bed on regular grids; subgrid cell-mean bed from sfincs_domain); on
+ ! subgrid models zb is overwritten with the water-level-dependent effective
+ ! bed every step (zb_effective), so the slopes must be taken now. The cap at +/- nonh_dzbmax
+ ! clips near-vertical wall steps that would otherwise produce an enormous
+ ! spurious w_b when a thin wall cell wets (2dt explicit instability); real
+ ! bed slopes are well below it.
+ !
+ allocate(nh_slbed(4, nrows))
+ nh_slbed = 0.0
+ !
+ block
+ integer :: iuv, nmn
+ real*4 :: slc
!
- ! Right
+ slc = nonh_dzbmax
+ if (slc <= 0.0) slc = huge(1.0) ! 0 = no cap
!
- if (nh_nm_index(2, irow) > 0) then
- !
- nmu = nm_index_of_row(nh_nm_index(2, irow))
- dzbdx(irow) = dzbdx(irow) + 0.5 * (zb(nmu) - zb(nm)) * dxrinv(1)
+ do irow = 1, nrows
!
- endif
- !
- ! Below
- !
- if (nh_nm_index(3, irow) > 0) then
+ nm = nm_index_of_row(irow)
!
- ndm = nm_index_of_row(nh_nm_index(3, irow))
- dzbdy(irow) = dzbdy(irow) + 0.5 * (zb(nm) - zb(ndm)) * dxrinv(1)
+ iuv = z_index_uv_md(nm)
+ if (iuv > 0) then
+ nmn = uv_index_z_nm(iuv)
+ nh_slbed(1, irow) = max(-slc, min(slc, (zb(nmn) - zb(nm)) * nh_dxr(irow)))
+ endif
!
- endif
- !
- ! Above
- !
- if (nh_nm_index(4, irow) > 0) then
+ iuv = z_index_uv_mu(nm)
+ if (iuv > 0) then
+ nmn = uv_index_z_nmu(iuv)
+ nh_slbed(2, irow) = max(-slc, min(slc, (zb(nm) - zb(nmn)) * nh_dxr(irow)))
+ endif
!
- num = nm_index_of_row(nh_nm_index(4, irow))
- dzbdy(irow) = dzbdy(irow) + 0.5 * (zb(num) - zb(nm)) * dxrinv(1)
+ if (nmax > 1) then
+ !
+ iuv = z_index_uv_nd(nm)
+ if (iuv > 0) then
+ nmn = uv_index_z_nm(iuv)
+ nh_slbed(3, irow) = max(-slc, min(slc, (zb(nmn) - zb(nm)) * nh_dyr(irow)))
+ endif
+ !
+ iuv = z_index_uv_nu(nm)
+ if (iuv > 0) then
+ nmn = uv_index_z_nmu(iuv)
+ nh_slbed(4, irow) = max(-slc, min(slc, (zb(nm) - zb(nmn)) * nh_dyr(irow)))
+ endif
+ !
+ endif
!
- endif
- !
- enddo
+ enddo
+ !
+ end block
+ !
+ ! Optional smoothing of the frozen bed slopes (nonh_slsmooth passes of a
+ ! [1 2 1]/4 stencil over the nonh neighbours, per direction). The single-layer
+ ! bottom kinematic condition w_b = u.d(zb)/dx makes the non-hydrostatic source
+ ! sensitive to the bed CURVATURE d2(zb)/dx2: at a slope break (e.g. the toe of
+ ! the conical island, where d(zb)/dx jumps over one cell) the per-cell w_b
+ ! jumps by ~ u.dx.d2(zb)/dx2, an impulsive grid-scale forcing that radiates a
+ ! trailing train of short dispersive waves. Smoothing ramps the slope over a
+ ! few cells, bounding the curvature (~ slope/(N.dx) instead of slope/dx) and
+ ! suppressing the wiggles at the source. Slopes are still frozen; this is the
+ ! curvature analogue of the nonh_dzbmax slope cap. Default 0 = no smoothing.
+ !
+ if (nonh_slsmooth > 0) then
+ block
+ real*4, dimension(:,:), allocatable :: sltmp
+ integer :: ps, ii, s, jw, je, js, jn
+ real*4 :: accs
+ allocate(sltmp(4, nrows))
+ do ps = 1, nonh_slsmooth
+ sltmp = nh_slbed
+ do ii = 1, nrows
+ jw = nh_nbw(ii) ; je = nh_nbe(ii)
+ do s = 1, 2 ! x-slopes smoothed along x-neighbours
+ accs = sltmp(s, ii)
+ if (jw > 0) accs = accs + 0.25 * (sltmp(s, jw) - sltmp(s, ii))
+ if (je > 0) accs = accs + 0.25 * (sltmp(s, je) - sltmp(s, ii))
+ nh_slbed(s, ii) = accs
+ enddo
+ js = nh_nbs(ii) ; jn = nh_nbn(ii)
+ do s = 3, 4 ! y-slopes smoothed along y-neighbours
+ accs = sltmp(s, ii)
+ if (js > 0) accs = accs + 0.25 * (sltmp(s, js) - sltmp(s, ii))
+ if (jn > 0) accs = accs + 0.25 * (sltmp(s, jn) - sltmp(s, ii))
+ nh_slbed(s, ii) = accs
+ enddo
+ enddo
+ enddo
+ deallocate(sltmp)
+ end block
+ endif
+ !
+ ! Open-boundary non-hydrostatic FADE-IN. Over nonh_fadein cells
+ ! inside each open (water level / outflow) boundary, ramp the non-hydrostatic
+ ! coupling from 0 (at the boundary) to full (nonh_fadein cells in). The incident
+ ! wave then enters hydrostatically at the boundary and the nonh turns on
+ ! gradually, avoiding the sharp pnh=0 -> full-nonh transition that reflects an
+ ! incoming dispersive wave. nh_fade is 1 everywhere by default (no fade).
+ !
+ if (nonh_fadein > 0) then
+ !
+ block
+ integer, dimension(:), allocatable :: idist
+ integer :: pass, nmL, nmR, rL, rR
+ !
+ allocate(idist(nrows))
+ idist = nonh_fadein + 1
+ do ip = 1, nhuv
+ inhuv = uv_index_of_nhuv(ip)
+ nmL = uv_index_z_nm(inhuv)
+ nmR = uv_index_z_nmu(inhuv)
+ rL = nh_faceL(ip)
+ rR = nh_faceR(ip)
+ if (rL > 0 .and. rR == 0) then
+ if (kcs(nmR) == 2 .or. kcs(nmR) == 3) idist(rL) = 0
+ endif
+ if (rR > 0 .and. rL == 0) then
+ if (kcs(nmL) == 2 .or. kcs(nmL) == 3) idist(rR) = 0
+ endif
+ enddo
+ do pass = 1, nonh_fadein
+ do ip = 1, nhuv
+ rL = nh_faceL(ip)
+ rR = nh_faceR(ip)
+ if (rL > 0 .and. rR > 0) then
+ if (idist(rL) + 1 < idist(rR)) idist(rR) = idist(rL) + 1
+ if (idist(rR) + 1 < idist(rL)) idist(rL) = idist(rR) + 1
+ endif
+ enddo
+ enddo
+ do irow = 1, nrows
+ nh_fade(irow) = min(1.0, real(idist(irow)) / real(nonh_fadein))
+ enddo
+ deallocate(idist)
+ end block
+ !
+ endif
+ !
+ ! Solver work arrays: assembled stencil, scaling and CG vectors (ghost
+ ! element 0). First-touch initialization (NUMA page placement matches the
+ ! static loop partition the per-step kernels use); ghost stays 0 forever.
+ !
+ allocate(nh_aw(nrows))
+ allocate(nh_ae(nrows))
+ allocate(nh_as(nrows))
+ allocate(nh_an(nrows))
+ allocate(nh_scl(nrows))
+ allocate(nh_sclo(0:nrows))
+ allocate(nh_y(0:nrows))
+ allocate(nh_r(0:nrows))
+ allocate(nh_p(0:nrows))
+ allocate(nh_s(0:nrows))
+ allocate(nh_w(0:nrows))
+ allocate(nh_b(0:nrows))
+ !
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_aw(i) = 0.0
+ nh_ae(i) = 0.0
+ nh_as(i) = 0.0
+ nh_an(i) = 0.0
+ nh_scl(i) = 1.0
+ nh_sclo(i) = 1.0
+ nh_y(i) = 0.0
+ nh_r(i) = 0.0
+ nh_p(i) = 0.0
+ nh_s(i) = 0.0
+ nh_w(i) = 0.0
+ nh_b(i) = 0.0
+ enddo
+ !$omp end parallel do
+ nh_sclo(0) = 0.0
+ nh_y(0) = 0.0
+ nh_r(0) = 0.0
+ nh_p(0) = 0.0
+ nh_s(0) = 0.0
+ nh_w(0) = 0.0
+ nh_b(0) = 0.0
+ !
+ ! Device residency for everything the per-step kernels read or write. The
+ ! global sfincs_data arrays (zs, zb, uv, q, kfuv, kcs, z_index_uv_*,
+ ! uv_index_z_*) are already on the device via initialize_openacc.
+ !
+ !$acc enter data copyin( nh_nbw, nh_nbe, nh_nbs, nh_nbn, &
+ !$acc nh_aw, nh_ae, nh_as, nh_an, nh_scl, nh_sclo, &
+ !$acc nh_y, nh_r, nh_p, nh_s, nh_w, nh_b, &
+ !$acc pnh, ws, wb, wb0, Dnm, nm_index_of_row, &
+ !$acc nh_dvert, nh_fade, nh_brfac, nh_bract, nh_dxr, nh_dyr, nh_slbed, &
+ !$acc nh_faceuv, nh_faceL, nh_faceR, nh_cellface, &
+ !$acc nh_cf, nh_cR, nh_cL, nh_hu )
+ !
+ write(logstr,'(a,i0,a,i0,a)')'Non-hydrostatic solver (assembled-stencil CG): ', nrows, ' cells, ', nhuv, ' faces'
+ call write_log(logstr, 0)
!
end subroutine
-
+
subroutine compute_nonhydrostatic(dt, tloop)
!
- ! Non-hydrostatic pressure correction on fluxes and velocities
+ ! Per-step non-hydrostatic pressure projection (see module header). Every
+ ! numbered kernel below is one !$omp / !$acc loop.
!
use sfincs_data
- use bicgstab_solver_ilu
!
implicit none
!
- integer :: count0
- integer :: count1
- integer :: count_rate
- integer :: count_max
+ integer :: count0, count1, count_rate, count_max
real :: tloop
- !
real*4 :: dt
!
- integer :: ip
- integer :: ipuv
- integer :: nm
- integer :: nmu
- integer :: nmd
- integer :: num
- integer :: ndm
- integer :: nmn
- integer :: j
- integer :: irow
- integer :: nhnm
- integer :: nhnmu
- integer :: iuv
- !
- integer :: iter
- !
- real*4 :: hu
- !
- real*4 :: Dnm1
- real*4 :: hnm
- real*4 :: Dnmu
- real*4 :: hnmu
- real*4 :: unh
- !
- real*4 :: hnb
- !
- real*4 :: dtover2rhodx2
- real*4 :: dtover2rhodx
- !
- real*4, dimension(npuv) :: AB
- real*4, dimension(:), allocatable :: QQ
- real*4, dimension(:), allocatable :: AA
- real*4 :: relres
+ integer :: i, ip, ipuv, nm, nmn, nmu, iuv, iuvl, iuvr, j, jl, jr, f, iter, ipass, cnt, ineighbour
+ real*4 :: dtrho, kbfac, hu, Dnm1, Dnmu, abf, fdep, tf, dval, sq, braw, accv, sumn, brrec, act
+ real*4 :: qr, ql, dzdt, wmax, breform, pcap, gf, pL, pR, unh
+ real*4 :: alpha, beta
+ real*8 :: gam8, del8, bn8, gamold8, alpha8, beta8, pap8, bnorm8
+ logical :: bndok
!
call system_clock(count0, count_rate, count_max)
!
- allocate(QQ(nrows))
- allocate(AA(nr_vals_in_matrix))
- !
- ! Compute AB (A and B)
+ if (nrows == 0) return
!
- AB = 0.0
- !
- !$omp parallel &
- !$omp private ( nm, irow )
- !$omp do schedule ( dynamic, 256 )
- do irow = 1, nrows
- !
- nm = nm_index_of_row(irow)
- !
- Dnm(irow) = max(zs(nm) - zb(nm), huthresh_nh)
- !
- enddo
- !$omp end do
- !$omp end parallel
- !
- !$omp parallel &
- !$omp private ( ip, ipuv, nm, nmu, num, hnm, hnmu, Dnm1, Dnmu )
- !$omp do schedule ( dynamic, 256 )
- do ip = 1, nhuv
- !
- ! Get water levels of neighboring cells
- !
- ! First get index of 'complete' uv array
- !
- ipuv = uv_index_of_nhuv(ip)
- !
- nm = uv_index_z_nm(ipuv)
- nmu = uv_index_z_nmu(ipuv)
- !
- if (kfuv(ipuv) == 1) then
- if (zs(nm) > zb(nm) + huthresh_nh .and. zs(nmu) > zb(nmu) + huthresh_nh) then
- !
- hnm = - zb(nm)
- Dnm1 = max(zs(nm) - zb(nm), huthresh_nh)
- !
- hnmu = - zb(nmu)
- Dnmu = max(zs(nmu) - zb(nmu), huthresh_nh)
- !
- AB(ip) = ( (zs(nmu) - hnmu) - (zs(nm) - hnm) ) / (Dnm1 + Dnmu)
- !
- endif
- endif
- !
- enddo
- !$omp end do
- !$omp end parallel
+ ! Keller-box vertical factor (nonh_disp): 2.0 = strict linear-pressure single layer;
+ ! lowering it (<2) strengthens the short-wave pressure feedback, flattening
+ ! c(k) toward Airy (~1.0 optimal). Default 1.0.
!
- ! Compute non-hydrostatic pressure by solving matrix AA * PP = QQ, where AA is a sparse matrix, PP is the nonh pressure, and QQ is the forcing
+ kbfac = nonh_disp
+ dtrho = dt / rhow
!
- AA = 0.0
- QQ = 0.0
+ ! breaking recovery per step: instant unless a reformation time is set
!
- dtover2rhodx2 = (dt * dxr2inv(1) / (2 * rhow))
+ brrec = 1.0
+ if (nonh_treform > 0.0) brrec = dt / nonh_treform
!
- ! Fill sparse matrix
+ ! 1) Layer depth per row and the vertical-acceleration diagonal.
!
- !$omp parallel &
- !$omp private ( ip, nm, j, nmd, nmu, ndm, num )
- !$omp do schedule ( dynamic, 256 )
- do irow = 1, nrows
- !
- ! Indices in nh uv array of neighboring uv points
- !
- nmd = nh_uv_index(1, irow)
- nmu = nh_uv_index(2, irow)
- ndm = nh_uv_index(3, irow)
- num = nh_uv_index(4, irow)
- !
- ! Left
- !
- if (nmd > 0) then
- !
- j = index_sparse_matrix(1, irow)
- !
- if (j>0) then
- !
- AA(j) = dtover2rhodx2 * (-1.0 + AB(nmd))
- !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, nm )
+ do i = 1, nrows
+ nm = nm_index_of_row(i)
+ Dnm(i) = max(zs(nm) - zb(nm), huthresh_nh)
+ nh_dvert(i) = kbfac * dt / (rhow * Dnm(i))
+ enddo
+ !$omp end parallel do
+ !
+ ! 1b) Breaking criterion -> nh_brfac (1 = full nonh, <1 = hydrostatic pnh=0).
+ ! Froude (NEOWAVE; Yamazaki, Kowalik & Cheung 2009, eqs 33/34) when
+ ! nonh_brfr > 0, else steepness/HFA (Smit, Zijlema & Stelling 2013) when
+ ! nonh_brsteep > 0.
+ !
+ if (nonh_brfr > 0.0) then
+ !
+ ! Froude criterion: a cell is set hydrostatic (pnh = 0) when the local flow
+ ! Froude number Fr = |U| / sqrt(g D) exceeds the onset nonh_brfr (~0.5), and
+ ! reactivated when Fr drops below 0.3*nonh_brfr -- NEOWAVE's hysteresis ratio.
+ ! Purely LOCAL (no smoothing / neighbour-spread / hole-filling needed) and
+ ! tracks a bore through its life.
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, nm, iuvl, iuvr, ql, qr, dzdt )
+ do i = 1, nrows
+ nm = nm_index_of_row(i)
+ ql = 0.0 ; qr = 0.0
+ iuvl = z_index_uv_md(nm) ; iuvr = z_index_uv_mu(nm)
+ if (iuvl > 0) ql = uv(iuvl)
+ if (iuvr > 0) qr = uv(iuvr)
+ dzdt = (0.5 * (ql + qr))**2 ! u_cell^2
+ if (nmax > 1) then
+ ql = 0.0 ; qr = 0.0
+ iuvl = z_index_uv_nd(nm) ; iuvr = z_index_uv_nu(nm)
+ if (iuvl > 0) ql = uv(iuvl)
+ if (iuvr > 0) qr = uv(iuvr)
+ dzdt = dzdt + (0.5 * (ql + qr))**2 ! + v_cell^2
endif
- !
- endif
- !
- ! Right
- !
- if (nmu > 0) then
- !
- j = index_sparse_matrix(2, irow)
- !
- if (j>0) then
- !
- AA(j) = dtover2rhodx2 * (-1.0 - AB(nmu))
- !
+ dzdt = sqrt(dzdt) / sqrt(g * Dnm(i)) ! Fr = |U| / sqrt(g D)
+ if (nh_brfac(i) >= 1.0) then ! was not breaking
+ if (dzdt > nonh_brfr) nh_brfac(i) = 0.0 ! onset
+ else ! was breaking -> release only when slow
+ if (dzdt < 0.3 * nonh_brfr) then
+ nh_brfac(i) = min(1.0, nh_brfac(i) + brrec) ! recover (gradually if nonh_treform > 0)
+ else
+ nh_brfac(i) = 0.0
+ endif
endif
- !
- endif
- !
- ! Bottom
- !
- if (ndm > 0) then
- !
- j = index_sparse_matrix(3, irow)
- !
- if (j>0) then
- !
- AA(j) = dtover2rhodx2 * (-1.0 + AB(ndm))
- !
+ enddo
+ !$omp end parallel do
+ !
+ elseif (nonh_brsteep > 0.0) then
+ !
+ ! Steepness/HFA criterion (gradual XBeach hydrostatic-front approximation).
+ ! Detect steepening from the surface rate of rise built from the flux
+ ! divergence: d(zs)/dx = dzdt / sqrt(g h), threshold nonh_brsteep. Only the
+ ! rising front (dzdt > 0) is reduced; a breaking cell stays breaking until
+ ! the surface falls.
+ !
+ breform = 0.25 * nonh_brsteep ! XBeach reformsteep default (neighbour-spread threshold)
+ !
+ ! (a) raw surface rate-of-rise dzdt = -div(q) into nh_w
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, nm, iuvl, iuvr, ql, qr, dzdt )
+ do i = 1, nrows
+ nm = nm_index_of_row(i)
+ iuvr = z_index_uv_mu(nm) ; iuvl = z_index_uv_md(nm)
+ qr = 0.0 ; ql = 0.0
+ if (iuvr > 0) qr = q(iuvr)
+ if (iuvl > 0) ql = q(iuvl)
+ dzdt = - (qr - ql) * nh_dxr(i)
+ if (nmax > 1) then
+ iuvr = z_index_uv_nu(nm) ; iuvl = z_index_uv_nd(nm)
+ qr = 0.0 ; ql = 0.0
+ if (iuvr > 0) qr = q(iuvr)
+ if (iuvl > 0) ql = q(iuvl)
+ dzdt = dzdt - (qr - ql) * nh_dyr(i)
endif
- !
- endif
- !
- ! Top
- !
- if (num > 0) then
- !
- j = index_sparse_matrix(4, irow)
- !
- if (j>0) then
- !
- AA(j) = dtover2rhodx2 * (-1.0 - AB(num))
- !
+ nh_w(i) = dzdt
+ enddo
+ !$omp end parallel do
+ !
+ ! (b) two [1 2 1]/4 smoothing passes (2dx noise -> connected front)
+ !
+ do ipass = 1, 2
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, accv, j )
+ do i = 1, nrows
+ accv = nh_w(i)
+ j = nh_nbw(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ j = nh_nbe(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ j = nh_nbs(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ j = nh_nbn(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ nh_p(i) = accv
+ enddo
+ !$omp end parallel do
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_w(i) = nh_p(i)
+ enddo
+ !$omp end parallel do
+ enddo
+ !
+ ! (c) snapshot last step's breaking state (no within-sweep contamination)
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_r(i) = nh_brfac(i)
+ enddo
+ !$omp end parallel do
+ !
+ ! (d) hysteretic state machine on the smoothed dzdt (onset / neighbour-spread /
+ ! release only when the surface falls)
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, dzdt, ineighbour, j, wmax )
+ do i = 1, nrows
+ dzdt = nh_w(i)
+ ineighbour = 0
+ j = nh_nbw(i) ; if (j > 0) then ; if (nh_r(j) < 1.0) ineighbour = 1 ; endif
+ j = nh_nbe(i) ; if (j > 0) then ; if (nh_r(j) < 1.0) ineighbour = 1 ; endif
+ j = nh_nbs(i) ; if (j > 0) then ; if (nh_r(j) < 1.0) ineighbour = 1 ; endif
+ j = nh_nbn(i) ; if (j > 0) then ; if (nh_r(j) < 1.0) ineighbour = 1 ; endif
+ wmax = sqrt(g * Dnm(i))
+ if (nh_r(i) >= 1.0) then ! was not breaking
+ if (dzdt > nonh_brsteep * wmax) then
+ nh_brfac(i) = 0.0
+ elseif (dzdt > breform * wmax .and. ineighbour == 1) then
+ nh_brfac(i) = 0.0
+ else
+ nh_brfac(i) = 1.0
+ endif
+ else ! was breaking
+ if (dzdt < 0.0) then
+ nh_brfac(i) = min(1.0, nh_r(i) + brrec) ! release only when surface falls (gradually if nonh_treform > 0)
+ else
+ nh_brfac(i) = 0.0
+ endif
+ endif
+ enddo
+ !$omp end parallel do
+ !
+ ! (e) close interior holes (a non-breaking cell flanked on both x-sides, or
+ ! both y-sides, by breaking cells is filled in; 3 passes)
+ !
+ do ipass = 1, 3
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_r(i) = nh_brfac(i)
+ enddo
+ !$omp end parallel do
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, ineighbour, jl, jr )
+ do i = 1, nrows
+ if (nh_r(i) >= 1.0) then
+ ineighbour = 0
+ jl = nh_nbw(i) ; if (jl > 0) then ; if (nh_r(jl) < 1.0) ineighbour = ineighbour + 1 ; endif
+ jr = nh_nbe(i) ; if (jr > 0) then ; if (nh_r(jr) < 1.0) ineighbour = ineighbour + 2 ; endif
+ if (ineighbour == 3) then
+ nh_brfac(i) = 0.0
+ endif
+ if (nmax > 1 .and. nh_brfac(i) >= 1.0) then
+ ineighbour = 0
+ jl = nh_nbs(i) ; if (jl > 0) then ; if (nh_r(jl) < 1.0) ineighbour = ineighbour + 1 ; endif
+ jr = nh_nbn(i) ; if (jr > 0) then ; if (nh_r(jr) < 1.0) ineighbour = ineighbour + 2 ; endif
+ if (ineighbour == 3) nh_brfac(i) = 0.0
+ endif
+ endif
+ enddo
+ !$omp end parallel do
+ enddo
+ !
+ else
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_brfac(i) = 1.0
+ enddo
+ !$omp end parallel do
+ !
+ endif
+ !
+ ! 1c) Breaking ACTIVITY for the pressure operator: nh_bract = nh_brfac,
+ ! optionally smoothed over the nonh neighbours (nonh_brsmooth passes of a
+ ! [1 2 1]/4 stencil). The binary flag switches the pressure off across a
+ ! single face, which prints a kink on the surface at the rear edge of the
+ ! breaking region; the smoothed activity ramps the pressure out over a few
+ ! cells instead. The hysteresis STATE stays in nh_brfac (unsmoothed).
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_bract(i) = nh_brfac(i)
+ enddo
+ !$omp end parallel do
+ !
+ do ipass = 1, nonh_brsmooth
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_w(i) = nh_bract(i)
+ enddo
+ !$omp end parallel do
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, accv, j )
+ do i = 1, nrows
+ accv = nh_w(i)
+ j = nh_nbw(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ j = nh_nbe(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ j = nh_nbs(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ j = nh_nbn(i) ; if (j > 0) accv = accv + 0.25 * (nh_w(j) - nh_w(i))
+ nh_bract(i) = accv
+ enddo
+ !$omp end parallel do
+ enddo
+ !
+ ! 2) Per-face gradient coefficients (wet/dry, layer term abf, conveyance depth,
+ ! open-boundary fade-in).
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( ip, ipuv, nm, nmu, hu, Dnm1, Dnmu, abf, fdep, bndok )
+ do ip = 1, nhuv
+ ipuv = nh_faceuv(ip)
+ nm = uv_index_z_nm(ipuv)
+ nmu = uv_index_z_nmu(ipuv)
+ nh_cR(ip) = 0.0
+ nh_cL(ip) = 0.0
+ nh_hu(ip) = 0.0
+ ! Open (kcs 2/3) boundary faces are skipped (handled hydrostatically);
+ ! closed walls / Neumann keep kfuv = 0 -> correct solid-wall BC.
+ bndok = (kcs(nm) == 1 .and. kcs(nmu) == 1)
+ if (kfuv(ipuv) == 1 .and. bndok) then
+ if (zs(nm) > zb(nm) + huthresh_nh .and. zs(nmu) > zb(nmu) + huthresh_nh) then
+ hu = max(zs(nm), zs(nmu)) - 0.5 * (zb(nm) + zb(nmu))
+ if (hu > huthresh_nh) then
+ Dnm1 = max(zs(nm) - zb(nm), huthresh_nh)
+ Dnmu = max(zs(nmu) - zb(nmu), huthresh_nh)
+ abf = ( (zs(nmu) + zb(nmu)) - (zs(nm) + zb(nm)) ) / (Dnm1 + Dnmu)
+ nh_cR(ip) = nh_cf(ip) * (1.0 + abf)
+ nh_cL(ip) = nh_cf(ip) * (abf - 1.0)
+ ! Conveyance depth q/uv keeps the corrected q and uv consistent
+ if (uv(ipuv) /= 0.0) then
+ nh_hu(ip) = q(ipuv) / uv(ipuv)
+ else
+ nh_hu(ip) = hu
+ endif
+ if (nonh_fadein > 0) then
+ fdep = 1.0
+ if (nh_faceL(ip) > 0) fdep = min(fdep, nh_fade(nh_faceL(ip)))
+ if (nh_faceR(ip) > 0) fdep = min(fdep, nh_fade(nh_faceR(ip)))
+ nh_cR(ip) = nh_cR(ip) * fdep
+ nh_cL(ip) = nh_cL(ip) * fdep
+ endif
+ endif
endif
- !
endif
- !
- ! Centre
- !
- j = index_sparse_matrix(5, irow)
- !
- ! if (j>0) then
- !
- AA(j) = 2 * dt / ( rhow * Dnm(irow)**2 )
- !
- ! endif
- !
- ! Forcing
- !
- QQ(irow) = 0.0
- !
- ! if (bnd(irow) == 0) then
- !
- QQ(irow) = - (ws(irow) + wb0(irow) - 2 * wb(irow)) / Dnm(irow)
- !
- if (nmd > 0) then
- !
- AA(j) = AA(j) + dtover2rhodx2 * (1.0 + AB(nmd))
- !
- QQ(irow) = QQ(irow) - (- uv(uv_index_of_nhuv(nmd))) * dxrinv(1)
- !
+ enddo
+ !$omp end parallel do
+ !
+ ! 3) ASSEMBLY (once per solve). Per row: raw off-diagonals a_ij = dtrho*hu*cR*cL
+ ! (one product per incident face, same both ways -> symmetric), the diagonal
+ ! dval = dvert + sum dtrho*hu*coef^2, the scaling scl = 1/sqrt(dval), the
+ ! breaking-masked scaling sclo, the SCALED right-hand side
+ ! b = scl * [ G^T (hu u*) - (ws + wb0 - 2 wb) ], and the warm start
+ ! y0 = pnh * sqrt(dval) (scaled previous pressure; 0 at breaking cells).
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, dval, braw, f, tf, sq, act )
+ do i = 1, nrows
+ dval = nh_dvert(i)
+ braw = - (ws(i) + wb0(i) - 2.0 * wb(i))
+ f = nh_cellface(1, i)
+ if (f > 0) then
+ tf = dtrho * nh_hu(f)
+ nh_aw(i) = tf * nh_cR(f) * nh_cL(f)
+ dval = dval + tf * nh_cR(f)**2
+ braw = braw + nh_cR(f) * nh_hu(f) * uv(nh_faceuv(f))
+ else
+ nh_aw(i) = 0.0
endif
- !
- if (nmu > 0) then
- !
- AA(j) = AA(j) + dtover2rhodx2 * (1.0 - AB(nmu))
- !
- QQ(irow) = QQ(irow) - (uv(uv_index_of_nhuv(nmu))) * dxrinv(1)
- !
+ f = nh_cellface(2, i)
+ if (f > 0) then
+ tf = dtrho * nh_hu(f)
+ nh_ae(i) = tf * nh_cL(f) * nh_cR(f)
+ dval = dval + tf * nh_cL(f)**2
+ braw = braw + nh_cL(f) * nh_hu(f) * uv(nh_faceuv(f))
+ else
+ nh_ae(i) = 0.0
endif
- !
- if (ndm > 0) then
- !
- AA(j) = AA(j) + dtover2rhodx2 * (1.0 + AB(ndm))
- !
- QQ(irow) = QQ(irow) - (- uv(uv_index_of_nhuv(ndm))) * dxrinv(1)
- !
+ f = nh_cellface(3, i)
+ if (f > 0) then
+ tf = dtrho * nh_hu(f)
+ nh_as(i) = tf * nh_cR(f) * nh_cL(f)
+ dval = dval + tf * nh_cR(f)**2
+ braw = braw + nh_cR(f) * nh_hu(f) * uv(nh_faceuv(f))
+ else
+ nh_as(i) = 0.0
endif
- !
- if (num > 0) then
- !
- AA(j) = AA(j) + dtover2rhodx2 * (1.0 - AB(num))
- !
- QQ(irow) = QQ(irow) - (uv(uv_index_of_nhuv(num))) * dxrinv(1)
- !
+ f = nh_cellface(4, i)
+ if (f > 0) then
+ tf = dtrho * nh_hu(f)
+ nh_an(i) = tf * nh_cL(f) * nh_cR(f)
+ dval = dval + tf * nh_cL(f)**2
+ braw = braw + nh_cL(f) * nh_hu(f) * uv(nh_faceuv(f))
+ else
+ nh_an(i) = 0.0
endif
- !
- ! endif
- !
+ sq = sqrt(dval) ! dval >= dvert > 0 always (Dnm capped)
+ nh_scl(i) = 1.0 / sq
+ !
+ ! Breaking activity act in [0,1]: act = 1 full nonh; act = 0 holds pnh = 0
+ ! (row+column zeroed via sclo, unit scaled diagonal). Fractional act (from
+ ! nonh_brsmooth / nonh_treform) ramps the pressure out smoothly; the
+ ! operator stays SPD (T A T + diag(1 - act^2), act <= 1).
+ !
+ act = nh_bract(i)
+ nh_sclo(i) = nh_scl(i) * act
+ nh_b(i) = braw * nh_scl(i) * act
+ nh_y(i) = pnh(i) * sq * act
enddo
- !$omp end do
- !$omp end parallel
- !
- ! Solve matrix
- !
- call bicgstab_solve(nrows, AA, col_idx, row_ptr, QQ, pnh, nh_tol, nh_itermax, iter, relres, .true.)
+ !$omp end parallel do
+ !
+ ! 3b) Scale the off-diagonals: a_ij <- a_ij * sclo(i) * sclo(j). Needs the
+ ! completed sclo of the NEIGHBOUR -> separate kernel. sclo(0) = 0 zeroes
+ ! coefficients to missing neighbours (the scaled diagonal is exactly 1).
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_aw(i) = nh_aw(i) * nh_sclo(i) * nh_sclo(nh_nbw(i))
+ nh_ae(i) = nh_ae(i) * nh_sclo(i) * nh_sclo(nh_nbe(i))
+ nh_as(i) = nh_as(i) * nh_sclo(i) * nh_sclo(nh_nbs(i))
+ nh_an(i) = nh_an(i) * nh_sclo(i) * nh_sclo(nh_nbn(i))
+ enddo
+ !$omp end parallel do
!
- ! Adjust fluxes
+ ! 4) CG solve of the scaled system (unit diagonal -> no preconditioner work).
+ ! Chronopoulos-Gear single-reduction CG: per iteration ONE fused
+ ! matvec+dots kernel and ONE fused vector-update kernel.
!
- dtover2rhodx = (dt * dxrinv(1) / (2 * rhow))
+ ! 4a) initial residual: w = A y0 (branch-free ghost matvec), r = b - w
!
- !$omp parallel &
- !$omp private ( ip, ipuv, nm, nmu, nhnm, nhnmu, hu, unh )
- !$omp do schedule ( dynamic, 256 )
- do ip = 1, nhuv
- !
- ipuv = uv_index_of_nhuv(ip)
- !
- if (kfuv(ipuv) == 1) then
- !
- ! Indices of neighbors in full zs array
- !
- nm = uv_index_z_nm(ipuv)
- nmu = uv_index_z_nmu(ipuv)
- !
- ! Indices of neighbors in nonh zs array
- !
- nhnm = row_index_of_nm(nm)
- nhnmu = row_index_of_nm(nmu)
- !
- hu = max(zs(nm), zs(nmu)) - 0.5 * (zb(nm) + zb(nmu))
- !
- if (hu > huthresh_nh) then
- !
- unh = - 1.0 * dtover2rhodx * ( AB(ip) * (pnh(nhnmu) + pnh(nhnm)) + pnh(nhnmu) - pnh(nhnm) )
- !
- ! Do some nudging to avoid 2dx waves
- !
- q(ipuv) = (1.0 - nh_fnudge) * q(ipuv) + nh_fnudge * (q(ipuv) + hu * unh)
- uv(ipuv) = (1.0 - nh_fnudge) * uv(ipuv) + nh_fnudge * (uv(ipuv) + unh)
- !
- endif
- !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_w(i) = nh_y(i) + nh_aw(i) * nh_y(nh_nbw(i)) &
+ + nh_ae(i) * nh_y(nh_nbe(i)) &
+ + nh_as(i) * nh_y(nh_nbs(i)) &
+ + nh_an(i) * nh_y(nh_nbn(i))
+ enddo
+ !$omp end parallel do
+ !
+ bn8 = 0.0d0
+ !$acc parallel loop default(present) reduction(+:bn8)
+ !$omp parallel do schedule ( static ) private ( i ) reduction(+:bn8)
+ do i = 1, nrows
+ nh_r(i) = nh_b(i) - nh_w(i)
+ nh_p(i) = 0.0
+ nh_s(i) = 0.0
+ bn8 = bn8 + nh_b(i) * nh_b(i)
+ enddo
+ !$omp end parallel do
+ !
+ bnorm8 = sqrt(bn8)
+ if (bnorm8 <= 0.0d0) bnorm8 = 1.0d0
+ !
+ ! 4b) iterate. Convergence on the scaled residual: ||S r|| / ||S b|| < nonh_tol.
+ !
+ iter = 0
+ gamold8 = 1.0d0
+ alpha8 = 1.0d0
+ !
+ do
+ !
+ ! kernel B: w = A r fused with gam = , del = (real*8 accumulators)
+ !
+ gam8 = 0.0d0
+ del8 = 0.0d0
+ !$acc parallel loop default(present) reduction(+:gam8,del8)
+ !$omp parallel do schedule ( static ) private ( i, accv ) reduction(+:gam8,del8)
+ do i = 1, nrows
+ accv = nh_r(i) + nh_aw(i) * nh_r(nh_nbw(i)) &
+ + nh_ae(i) * nh_r(nh_nbe(i)) &
+ + nh_as(i) * nh_r(nh_nbs(i)) &
+ + nh_an(i) * nh_r(nh_nbn(i))
+ nh_w(i) = accv
+ gam8 = gam8 + nh_r(i) * nh_r(i)
+ del8 = del8 + nh_r(i) * accv
+ enddo
+ !$omp end parallel do
+ !
+ if (sqrt(gam8) < nonh_tol * bnorm8) exit
+ if (iter >= nonh_itermax) exit
+ !
+ ! scalar recurrences (Chronopoulos-Gear): beta, pAp, alpha
+ !
+ if (iter == 0) then
+ beta8 = 0.0d0
+ pap8 = del8
+ else
+ beta8 = gam8 / gamold8
+ pap8 = del8 - beta8 * gam8 / alpha8 ! =
endif
+ if (pap8 <= 0.0d0) exit ! breakdown guard (SPD -> pAp > 0)
+ alpha8 = gam8 / pap8
+ gamold8 = gam8
+ alpha = real(alpha8, 4)
+ beta = real(beta8, 4)
+ !
+ ! kernel A: p = r + beta p ; s = w + beta s ; y += alpha p ; r -= alpha s
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_p(i) = nh_r(i) + beta * nh_p(i)
+ nh_s(i) = nh_w(i) + beta * nh_s(i)
+ nh_y(i) = nh_y(i) + alpha * nh_p(i)
+ nh_r(i) = nh_r(i) - alpha * nh_s(i)
+ enddo
+ !$omp end parallel do
+ !
+ iter = iter + 1
!
- enddo
- !$omp end do
- !$omp end parallel
+ enddo
!
- ! Update vertical velocity ws and wb
+ ! diagnostics (reported in the end-of-run timing summary)
!
- !$omp parallel &
- !$omp private ( nm, iuv, irow, nmn, hnm, hnb)
- !$omp do schedule ( dynamic, 256 )
- do irow = 1, nrows
- !
- ! Copy wb0 from previous time step
- !
- wb0(irow) = wb(irow)
- !
- wb(irow) = 0.0
- !
- nm = nm_index_of_row(irow)
- hnm = zs(nm) - zb(nm)
- !
- ! This will not yet work for quadtree !
- !
- ! Indices of neighboring cells
- !
- ! Left
- !
- iuv = z_index_uv_md(nm) ! uv index
- !
- if (kfuv(iuv) > 0) then
- nmn = uv_index_z_nm(iuv) ! nm index of neighbor
- hnb = zs(nmn) - zb(nmn) ! water depth at neighbor
- wb(irow) = wb(irow) - 0.5 * uv(iuv) * (hnm - hnb) * dxrinv(1)
+ nh_iter_total = nh_iter_total + iter
+ nh_solve_count = nh_solve_count + 1
+ if (iter > nh_iter_max) nh_iter_max = iter
+ !
+ ! 4c) unscale: pnh = S y (breaking rows have y = 0 -> pnh = 0 held exactly)
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ pnh(i) = nh_scl(i) * nh_y(i)
+ enddo
+ !$omp end parallel do
+ !
+ ! 5b) Optional global 2dx filter on pnh (nonh_filter > 0): one Jacobi
+ ! neighbour-mean pass; snapshot into nh_w first (race-free). Damps the
+ ! persistent grid-scale mode without touching the resolved smooth pressure.
+ !
+ if (nonh_filter > 0.0) then
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_w(i) = pnh(i)
+ enddo
+ !$omp end parallel do
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, sumn, cnt, j )
+ do i = 1, nrows
+ sumn = 0.0 ; cnt = 0
+ j = nh_nbw(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ j = nh_nbe(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ j = nh_nbs(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ j = nh_nbn(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ if (cnt > 0) pnh(i) = (1.0 - nonh_filter) * nh_w(i) + nonh_filter * sumn / real(cnt)
+ enddo
+ !$omp end parallel do
+ endif
+ !
+ ! 5c) Localized 2dx smoothing in the marginal-nonh zones (fade-in zone and/or
+ ! shallow run-up water, D < nonh_smoothdep); zero in the resolved interior.
+ !
+ if (nonh_smoothbnd > 0.0 .and. (nonh_fadein > 0 .or. nonh_smoothdep > 0.0)) then
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i )
+ do i = 1, nrows
+ nh_w(i) = pnh(i)
+ enddo
+ !$omp end parallel do
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, sumn, cnt, j, fdep, gf )
+ do i = 1, nrows
+ sumn = 0.0 ; cnt = 0
+ j = nh_nbw(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ j = nh_nbe(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ j = nh_nbs(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ j = nh_nbn(i) ; if (j > 0) then ; sumn = sumn + nh_w(j) ; cnt = cnt + 1 ; endif
+ if (cnt > 0) then
+ fdep = 0.0
+ if (nonh_fadein > 0) fdep = 1.0 - nh_fade(i)
+ if (nonh_smoothdep > 0.0) fdep = max(fdep, (nonh_smoothdep - Dnm(i)) / nonh_smoothdep)
+ gf = nonh_smoothbnd * max(0.0, min(1.0, fdep))
+ pnh(i) = (1.0 - gf) * nh_w(i) + gf * sumn / real(cnt)
+ endif
+ enddo
+ !$omp end parallel do
+ endif
+ !
+ ! 5e) Depth limiter: |pnh| <= nonh_pmax * rho g H. The total bed pressure
+ ! rho*g*H + pnh then stays >= 0 for nonh_pmax <= 1 (no "suction"); resolved
+ ! waves have |pnh| << rho*g*H so are untouched -- this only clips spikes
+ ! at steep wall run-up, breaking fronts and grid-scale modes.
+ !
+ if (nonh_pmax > 0.0) then
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, pcap )
+ do i = 1, nrows
+ pcap = nonh_pmax * rhow * g * Dnm(i)
+ pnh(i) = max(-pcap, min(pcap, pnh(i)))
+ enddo
+ !$omp end parallel do
+ endif
+ !
+ ! 6) Correct fluxes / velocities with the pressure gradient -(dt/rho) G pnh.
+ ! Each face writes only its own uv/q point -> race-free.
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( ip, ipuv, pL, pR, gf, unh )
+ do ip = 1, nhuv
+ if (nh_cR(ip) == 0.0 .and. nh_cL(ip) == 0.0) cycle
+ ipuv = nh_faceuv(ip)
+ pR = 0.0
+ pL = 0.0
+ if (nh_faceR(ip) > 0) pR = pnh(nh_faceR(ip))
+ if (nh_faceL(ip) > 0) pL = pnh(nh_faceL(ip))
+ gf = nh_cR(ip) * pR + nh_cL(ip) * pL
+ unh = - dtrho * gf
+ uv(ipuv) = (1.0 - nonh_fnudge) * uv(ipuv) + nonh_fnudge * (uv(ipuv) + unh)
+ q(ipuv) = (1.0 - nonh_fnudge) * q(ipuv) + nonh_fnudge * (q(ipuv) + nh_hu(ip) * unh)
+ enddo
+ !$omp end parallel do
+ !
+ ! 7) Update surface/bottom vertical velocities for the next step's forcing
+ ! (bottom kinematic condition w_b = u . d(zb)/dx with the slopes FROZEN at
+ ! initialization in nh_slbed; the wet checks use the current, possibly
+ ! effective, zb -- they are depth checks).
+ !
+ !$acc parallel loop default(present)
+ !$omp parallel do schedule ( static ) private ( i, nm, iuv, nmn )
+ do i = 1, nrows
+ nm = nm_index_of_row(i)
+ !
+ ! Non-hydrostatically dry cell: carry no nonh state (clean re-wetting).
+ !
+ if (zs(nm) - zb(nm) <= huthresh_nh) then
+ wb0(i) = 0.0
+ wb(i) = 0.0
+ ws(i) = 0.0
+ pnh(i) = 0.0
+ nh_brfac(i) = 1.0
+ cycle
endif
!
- ! Right
+ wb0(i) = wb(i)
+ wb(i) = 0.0
!
- iuv = z_index_uv_mu(nm) ! uv index
+ ! w_b over faces to non-hydrostatically wet neighbours only.
!
+ iuv = z_index_uv_md(nm)
if (kfuv(iuv) > 0) then
- nmn = uv_index_z_nmu(iuv) ! nm index of neighbor
- hnb = zs(nmn) - zb(nmn) ! water depth at neighbor
- wb(irow) = wb(irow) - 0.5 * uv(iuv) * (hnb - hnm) * dxrinv(1)
+ nmn = uv_index_z_nm(iuv)
+ if (zs(nmn) - zb(nmn) > huthresh_nh) then
+ wb(i) = wb(i) - 0.5 * uv(iuv) * nh_slbed(1, i)
+ endif
endif
- !
- ! Bottom
- !
- iuv = z_index_uv_nd(nm) ! uv index
- !
+ iuv = z_index_uv_mu(nm)
if (kfuv(iuv) > 0) then
- nmn = uv_index_z_nm(iuv) ! nm index of neighbor
- hnb = zs(nmn) - zb(nmn) ! water depth at neighbor
- wb(irow) = wb(irow) - 0.5 * uv(iuv) * (hnm - hnb) * dyrinv(1)
+ nmn = uv_index_z_nmu(iuv)
+ if (zs(nmn) - zb(nmn) > huthresh_nh) then
+ wb(i) = wb(i) - 0.5 * uv(iuv) * nh_slbed(2, i)
+ endif
endif
- !
- ! Top
- !
- iuv = z_index_uv_nu(nm) ! uv index
- !
- if (kfuv(iuv) > 0) then
- nmn = uv_index_z_nmu(iuv) ! nm index of neighbor
- hnb = zs(nmn) - zb(nmn) ! water depth at neighbor
- wb(irow) = wb(irow) - 0.5 * uv(iuv) * (hnb - hnm) * dyrinv(1)
+ if (nmax > 1) then
+ iuv = z_index_uv_nd(nm)
+ if (kfuv(iuv) > 0) then
+ nmn = uv_index_z_nm(iuv)
+ if (zs(nmn) - zb(nmn) > huthresh_nh) then
+ wb(i) = wb(i) - 0.5 * uv(iuv) * nh_slbed(3, i)
+ endif
+ endif
+ iuv = z_index_uv_nu(nm)
+ if (kfuv(iuv) > 0) then
+ nmn = uv_index_z_nmu(iuv)
+ if (zs(nmn) - zb(nmn) > huthresh_nh) then
+ wb(i) = wb(i) - 0.5 * uv(iuv) * nh_slbed(4, i)
+ endif
+ endif
endif
!
- ws(irow) = ws(irow) - (wb(irow) - wb0(irow)) + (2 * dt / (rhow * Dnm(irow))) * pnh(irow) ! this is ws m+1 in the next time step
- !
- enddo
- !$omp end do
- !$omp end parallel
+ ws(i) = ws(i) - (wb(i) - wb0(i)) + (kbfac * dt / (rhow * Dnm(i))) * pnh(i)
+ enddo
+ !$omp end parallel do
+ !
+ ! Optional moving-bed source: add d(zb)/dt = dzbext/dt to the bottom kinematic
+ ! w_b (full condition w_b = d(zb)/dt + u.d(zb)/dx), so a moving seafloor
+ ! radiates a depth-filtered (Kajiura-like) surface response instead of being
+ ! stamped onto zs. Kept in a SEPARATE flag-guarded loop so the step-7 loop
+ ! above is byte-for-byte the original when the feature is off. Adding the term
+ ! here (wb += d(zb)/dt, ws -= d(zb)/dt) is algebraically identical to carrying
+ ! it inside step 7 (it cancels through the ws = ws - (wb - wb0) update).
+ if (nonh_movingbed .and. use_dzbext .and. dt > 0.0) then
+ !$omp parallel do schedule ( static ) private ( i, nm )
+ do i = 1, nrows
+ nm = nm_index_of_row(i)
+ if (zs(nm) - zb(nm) <= huthresh_nh) cycle
+ wb(i) = wb(i) + dzbext(nm) / dt
+ ws(i) = ws(i) - dzbext(nm) / dt
+ enddo
+ !$omp end parallel do
+ endif
!
call system_clock(count1, count_rate, count_max)
tloop = tloop + 1.0*(count1 - count0)/count_rate
!
- end subroutine
+ end subroutine
end module
diff --git a/source/src/sfincs_output.f90 b/source/src/sfincs_output.f90
index bacd6750f..77a78dbf4 100644
--- a/source/src/sfincs_output.f90
+++ b/source/src/sfincs_output.f90
@@ -80,6 +80,7 @@ subroutine initialize_output(tmapout,tmaxout,thisout, trstout)
subroutine write_output(t,write_map,write_his,write_max,write_rst,ntmapout,ntmaxout,nthisout,tloop)
!
use sfincs_data
+ use sfincs_nonhydrostatic, only: pnh
!
implicit none
!
@@ -110,6 +111,15 @@ subroutine write_output(t,write_map,write_his,write_max,write_rst,ntmapout,ntmax
!
!$acc update host(zs)
!
+ if (nonhydrostatic .and. write_map) then
+ !
+ ! Non-hydrostatic pressure is written to the map file (device-resident
+ ! in GPU builds; comment / no-op in CPU builds)
+ !
+ !$acc update host(pnh) if_present
+ !
+ endif
+ !
if (store_cumulative_precipitation) then
!
!$acc update host(prcp)
diff --git a/source/third_party_open/netcdf/netcdf-fortran-4.6.1/netcdff_c.vcxproj b/source/third_party_open/netcdf/netcdf-fortran-4.6.1/netcdff_c.vcxproj
index bb14a4db0..e79c9c068 100644
--- a/source/third_party_open/netcdf/netcdf-fortran-4.6.1/netcdff_c.vcxproj
+++ b/source/third_party_open/netcdf/netcdf-fortran-4.6.1/netcdff_c.vcxproj
@@ -25,12 +25,12 @@
StaticLibrary
MultiByte
- v143
+ v145
StaticLibrary
MultiByte
- v143
+ v145