From a6f3ee92a9c838144471aeaf49359b7ea65e1f6a Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 13 Feb 2026 11:16:57 -0700 Subject: [PATCH 01/12] Normals are calculated together with the gradients --- physics/level_set_mod.F90 | 74 +++++++++++++++++++++++++++++---------- 1 file changed, 56 insertions(+), 18 deletions(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 0f3d35e..c5e96c4 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -986,7 +986,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm real, parameter :: EPS = epsilon (0.0), TOL = 100.0 * EPS real :: difflx, diffly, diffrx, diffry, diffcx, diffcy, & - diff2x, diff2y, grad, & + diff2x, diff2y, grad, gradx, grady, & scale, nvx, nvy, a_valor, signo_x, signo_y, threshold_hll, & threshold_hlu, threshold_av, fire_viscosity_var integer :: i, j @@ -1024,40 +1024,67 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) else select case (fire_upwinding) - ! none case (0) + ! none grad = sqrt (diffcx ** 2 + diffcy ** 2) - ! standard + scale = sqrt (grad ** 2.0 + EPS) + nvx = diffcx / scale + nvy = diffcy / scale + case (1) + ! standard diff2x = Select_upwind (difflx, diffrx) diff2y = Select_upwind (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! godunov per osher/fedkiw + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case (2) + ! godunov per osher/fedkiw diff2x = Select_godunov (difflx, diffrx) diff2y = Select_godunov (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! ENO1 + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case (3) + ! ENO1 diff2x = Select_eno (difflx, diffrx) diff2y = Select_eno (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! Sethian - twice stronger pushdown of bumps + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case(4) - grad = sqrt (max (difflx, 0.0) ** 2 + min (diffrx, 0.0) ** 2 & - + max (diffly, 0.0) ** 2 + min(diffry, 0.0) ** 2) - ! 2nd order + ! Sethian - twice stronger pushdown of bumps + grad = sqrt (max (difflx, 0.0) ** 2 + min (diffrx, 0.0) ** 2 + & + max (diffly, 0.0) ** 2 + min(diffry, 0.0) ** 2) + + scale = sqrt (grad ** 2.0 + EPS) + gradx = max(difflx,0.0) - min(diffrx,0.0) + grady = max(diffly,0.0) - min(diffry,0.0) + nvx = gradx / scale + nvy = grady / scale + case(5) + ! 2nd order diff2x = Select_2nd (dx, lfn(i, j), lfn(i - 1, j), lfn(i + 1, j)) diff2y = Select_2nd (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j + 1)) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! WENO3 + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case(6) + ! WENO3 a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), & @@ -1070,8 +1097,12 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm lfn(i, j + 1), lfn(i, j + 2), signo_y) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! WENO5 + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case(7) + ! WENO5 a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j)+ & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) @@ -1082,8 +1113,12 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm lfn(i ,j - 3), lfn(i, j + 1), lfn(i, j + 2), lfn(i, j + 3), signo_y) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - ! WENO3/ENO1 + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case(8) + ! WENO3/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) @@ -1102,8 +1137,12 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) end if - ! WENO5/ENO1 + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case(9) + ! WENO5/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & Select_4th (dy,lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) @@ -1122,6 +1161,10 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) end if + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case default !$omp critical write (msg, '(a, i2)') 'Unknown upwinding option in level set : ', fire_upwinding @@ -1131,11 +1174,6 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm end select end if - ! Calc normal - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - ! Get rate of spread from wind speed and slope ros(i, j) = ros_model%Calc_ros (ifms, ifme, jfms, jfme, i, j, & nvx, nvy, uf(i, j), vf(i, j), dzdxf(i, j), dzdyf(i, j)) From 548e2d1b58d8cec8521b4a3615e094839ccb17a7 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 13 Feb 2026 12:08:32 -0700 Subject: [PATCH 02/12] Adding upwinding with transition xone from weno to eno --- physics/level_set_mod.F90 | 45 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index c5e96c4..563d6a5 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -18,6 +18,7 @@ module level_set_mod use state_mod, only: state_fire_t, N_POINTS_IN_HALO use ignition_line_mod, only : ignition_line_t use ros_mod, only : ros_t + use constants_mod, only : PI #ifdef DM_PARALLEL use mpi_mod, only : Do_halo_exchange, Sum_across_mpi_tasks, Max_across_mpi_tasks @@ -986,7 +987,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm real, parameter :: EPS = epsilon (0.0), TOL = 100.0 * EPS real :: difflx, diffly, diffrx, diffry, diffcx, diffcy, & - diff2x, diff2y, grad, gradx, grady, & + diff2x, diff2y, grad, gradx, grady, mask, diff2x_eno, diff2y_eno, & + diff2x_weno, diff2y_weno, transition_width, band_width, & scale, nvx, nvy, a_valor, signo_x, signo_y, threshold_hll, & threshold_hlu, threshold_av, fire_viscosity_var integer :: i, j @@ -1165,6 +1167,47 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvx = diff2x / scale nvy = diff2y / scale + case(10) + ! WENO5/ENO1 + blending zone + band_width = 2 * threshold_hlu + transition_width = threshold_hlu + if (abs(lfn(i, j)) <= band_width) then + if (abs(lfn(i, j)) <= band_width - transition_width) then + mask = 1.0 + else + mask = 0.5 * (1.0 + cos (PI * (abs (lfn(i, j)) -(band_width - transition_width)) / transition_width)) + end if + else + mask = 0.0 + end if + + diff2x_eno = Select_eno (difflx, diffrx) + diff2y_eno = Select_eno (diffly, diffry) + + if (mask > 0.0) then + a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & + Select_4th (dy,lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) + signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), & + lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) + signo_y = a_valor * Select_4th (dy, lfn(i, j), lfn(i, j - 1), & + lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) + diff2x_weno = Select_weno5 (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), & + lfn(i - 3, j), lfn(i + 1, j), lfn(i + 2, j), lfn(i + 3, j), signo_x) + diff2y_weno = Select_weno5 (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), & + lfn(i, j - 3), lfn(i, j + 1), lfn(i, j + 2), lfn(i, j + 3), signo_y) + diff2x = mask * diff2x_weno + (1.0 - mask) * diff2x_eno + diff2y = mask * diff2y_weno + (1.0 - mask) * diff2y_eno + else + diff2x = diff2x_eno + diff2y = diff2y_eno + end if + + grad = sqrt (diff2x * diff2x + diff2y * diff2y) + + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + case default !$omp critical write (msg, '(a, i2)') 'Unknown upwinding option in level set : ', fire_upwinding From 1ea8d536d3ab887aab30736a3901da81489b7ec0 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 13 Feb 2026 16:33:11 -0700 Subject: [PATCH 03/12] Add missing factor in central differences calculation --- physics/level_set_mod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 563d6a5..df188fb 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -1015,9 +1015,9 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm difflx = (lfn(i, j) - lfn(i - 1, j)) / dx diffry = (lfn(i, j + 1) - lfn(i, j)) / dy diffly = (lfn(i, j) - lfn(i, j - 1)) / dy - ! twice central difference - diffcx = difflx + diffrx - diffcy = diffly + diffry + ! central difference + diffcx = 0.5 * (difflx + diffrx) + diffcy = 0.5 * (diffly + diffry) ! use eno1 near domain boundaries if (i < ids + BDY_ENO1 .or. i > ide - BDY_ENO1 .or. & j < jds + BDY_ENO1 .or. j > jde - BDY_ENO1) then From 0fb2d7641cd0653cdc82f12b161fb5ca30694d58 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 13 Feb 2026 16:41:59 -0700 Subject: [PATCH 04/12] fixes calculation of upwinding 4 --- physics/level_set_mod.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index df188fb..531d514 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -987,7 +987,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm real, parameter :: EPS = epsilon (0.0), TOL = 100.0 * EPS real :: difflx, diffly, diffrx, diffry, diffcx, diffcy, & - diff2x, diff2y, grad, gradx, grady, mask, diff2x_eno, diff2y_eno, & + diff2x, diff2y, grad, mask, diff2x_eno, diff2y_eno, & diff2x_weno, diff2y_weno, transition_width, band_width, & scale, nvx, nvy, a_valor, signo_x, signo_y, threshold_hll, & threshold_hlu, threshold_av, fire_viscosity_var @@ -1067,13 +1067,13 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm case(4) ! Sethian - twice stronger pushdown of bumps grad = sqrt (max (difflx, 0.0) ** 2 + min (diffrx, 0.0) ** 2 + & - max (diffly, 0.0) ** 2 + min(diffry, 0.0) ** 2) + max (diffly, 0.0) ** 2 + min (diffry, 0.0) ** 2) scale = sqrt (grad ** 2.0 + EPS) - gradx = max(difflx,0.0) - min(diffrx,0.0) - grady = max(diffly,0.0) - min(diffry,0.0) - nvx = gradx / scale - nvy = grady / scale + diff2x = max (difflx, 0.0) - min (diffrx, 0.0) + diff2y = max (diffly, 0.0) - min (diffry, 0.0) + nvx = diff2x / scale + nvy = diff2y / scale case(5) ! 2nd order From 0ff3803b1609fe106bcdb067ae736fa5a27ce4d9 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 13 Feb 2026 16:54:35 -0700 Subject: [PATCH 05/12] Exposes the gradient to CFL condition for no upwinding case --- physics/level_set_mod.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 531d514..4ec64a9 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -1028,11 +1028,13 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm select case (fire_upwinding) case (0) ! none - grad = sqrt (diffcx ** 2 + diffcy ** 2) + diff2x = diffcx + diff2y = diffcy + grad = sqrt (diff2x * diff2x + diff2y * diff2y) scale = sqrt (grad ** 2.0 + EPS) - nvx = diffcx / scale - nvy = diffcy / scale + nvx = diff2x / scale + nvy = diff2y / scale case (1) ! standard From 8300b832144efdf2eb948fb684b41700001af5e7 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 13 Feb 2026 17:01:12 -0700 Subject: [PATCH 06/12] Calc normals near the boundary too --- physics/level_set_mod.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 4ec64a9..f76e561 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -1015,15 +1015,21 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm difflx = (lfn(i, j) - lfn(i - 1, j)) / dx diffry = (lfn(i, j + 1) - lfn(i, j)) / dy diffly = (lfn(i, j) - lfn(i, j - 1)) / dy - ! central difference + + ! central differences diffcx = 0.5 * (difflx + diffrx) diffcy = 0.5 * (diffly + diffry) + ! use eno1 near domain boundaries if (i < ids + BDY_ENO1 .or. i > ide - BDY_ENO1 .or. & j < jds + BDY_ENO1 .or. j > jde - BDY_ENO1) then diff2x = Select_eno (difflx, diffrx) diff2y = Select_eno (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) + + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale else select case (fire_upwinding) case (0) From 52f170c351c8331d79c2afb50ca59f990dc34783 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Sat, 14 Feb 2026 07:17:29 -0700 Subject: [PATCH 07/12] Impro comments --- physics/level_set_mod.F90 | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index f76e561..0339454 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -1043,7 +1043,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvy = diff2y / scale case (1) - ! standard + ! First-order sign-consistent upwind scheme (component wise). Standard diff2x = Select_upwind (difflx, diffrx) diff2y = Select_upwind (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) @@ -1053,7 +1053,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvy = diff2y / scale case (2) - ! godunov per osher/fedkiw + ! Component-wise Godunov upwind (1st order) derivative selection diff2x = Select_godunov (difflx, diffrx) diff2y = Select_godunov (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) @@ -1063,7 +1063,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvy = diff2y / scale case (3) - ! ENO1 + ! First order ENO (minmod limited) componenet-wise gradient diff2x = Select_eno (difflx, diffrx) diff2y = Select_eno (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) @@ -1072,19 +1072,19 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvx = diff2x / scale nvy = diff2y / scale - case(4) - ! Sethian - twice stronger pushdown of bumps + case (4) + ! Classical Godunov discretization of the Hamiltonian grad = sqrt (max (difflx, 0.0) ** 2 + min (diffrx, 0.0) ** 2 + & max (diffly, 0.0) ** 2 + min (diffry, 0.0) ** 2) - - scale = sqrt (grad ** 2.0 + EPS) diff2x = max (difflx, 0.0) - min (diffrx, 0.0) diff2y = max (diffly, 0.0) - min (diffry, 0.0) + + scale = sqrt (grad ** 2.0 + EPS) nvx = diff2x / scale nvy = diff2y / scale - case(5) - ! 2nd order + case (5) + ! 2nd order central differences diff2x = Select_2nd (dx, lfn(i, j), lfn(i - 1, j), lfn(i + 1, j)) diff2y = Select_2nd (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j + 1)) grad = sqrt (diff2x * diff2x + diff2y * diff2y) @@ -1093,8 +1093,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvx = diff2x / scale nvy = diff2y / scale - case(6) - ! WENO3 + case (6) + ! 3rd order WENO a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), & @@ -1111,8 +1111,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvx = diff2x / scale nvy = diff2y / scale - case(7) - ! WENO5 + case (7) + ! 5th order WENO a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j)+ & Select_4th (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j - 2), lfn(i, j + 1), lfn(i, j + 2)) * vf(i, j) signo_x = a_valor * Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) @@ -1127,7 +1127,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvx = diff2x / scale nvy = diff2y / scale - case(8) + case (8) ! WENO3/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & @@ -1151,7 +1151,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvx = diff2x / scale nvy = diff2y / scale - case(9) + case (9) ! WENO5/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & @@ -1175,7 +1175,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm nvx = diff2x / scale nvy = diff2y / scale - case(10) + case (10) ! WENO5/ENO1 + blending zone band_width = 2 * threshold_hlu transition_width = threshold_hlu From d47373b83e38ebd8263c6e0d73d75715f3d6adf0 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Sat, 14 Feb 2026 08:19:15 -0700 Subject: [PATCH 08/12] Normals calculated outside of the upwinding selection --- physics/level_set_mod.F90 | 52 ++++----------------------------------- 1 file changed, 5 insertions(+), 47 deletions(-) diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 0339454..94809e2 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -1027,9 +1027,6 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm diff2y = Select_eno (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale else select case (fire_upwinding) case (0) @@ -1038,40 +1035,24 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm diff2y = diffcy grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (1) ! First-order sign-consistent upwind scheme (component wise). Standard diff2x = Select_upwind (difflx, diffrx) diff2y = Select_upwind (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (2) ! Component-wise Godunov upwind (1st order) derivative selection diff2x = Select_godunov (difflx, diffrx) diff2y = Select_godunov (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (3) ! First order ENO (minmod limited) componenet-wise gradient diff2x = Select_eno (difflx, diffrx) diff2y = Select_eno (diffly, diffry) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (4) ! Classical Godunov discretization of the Hamiltonian grad = sqrt (max (difflx, 0.0) ** 2 + min (diffrx, 0.0) ** 2 + & @@ -1079,20 +1060,12 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm diff2x = max (difflx, 0.0) - min (diffrx, 0.0) diff2y = max (diffly, 0.0) - min (diffry, 0.0) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (5) ! 2nd order central differences diff2x = Select_2nd (dx, lfn(i, j), lfn(i - 1, j), lfn(i + 1, j)) diff2y = Select_2nd (dy, lfn(i, j), lfn(i, j - 1), lfn(i, j + 1)) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (6) ! 3rd order WENO a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j) + & @@ -1107,10 +1080,6 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm lfn(i, j + 1), lfn(i, j + 2), signo_y) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (7) ! 5th order WENO a_valor = Select_4th (dx, lfn(i, j), lfn(i - 1, j), lfn(i - 2, j), lfn(i + 1, j), lfn(i + 2, j)) * uf(i, j)+ & @@ -1123,10 +1092,6 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm lfn(i ,j - 3), lfn(i, j + 1), lfn(i, j + 2), lfn(i, j + 3), signo_y) grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (8) ! WENO3/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then @@ -1147,10 +1112,6 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) end if - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (9) ! WENO5/ENO1 if (abs (lfn(i, j)) < threshold_hlu) then @@ -1171,10 +1132,6 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) end if - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case (10) ! WENO5/ENO1 + blending zone band_width = 2 * threshold_hlu @@ -1212,10 +1169,6 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm grad = sqrt (diff2x * diff2x + diff2y * diff2y) - scale = sqrt (grad ** 2.0 + EPS) - nvx = diff2x / scale - nvy = diff2y / scale - case default !$omp critical write (msg, '(a, i2)') 'Unknown upwinding option in level set : ', fire_upwinding @@ -1225,6 +1178,11 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm end select end if + ! Calc normal + scale = sqrt (grad ** 2.0 + EPS) + nvx = diff2x / scale + nvy = diff2y / scale + ! Get rate of spread from wind speed and slope ros(i, j) = ros_model%Calc_ros (ifms, ifme, jfms, jfme, i, j, & nvx, nvy, uf(i, j), vf(i, j), dzdxf(i, j), dzdyf(i, j)) From 79883686a3a579f82d97841c449eb9dbf7a73baf Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Sun, 15 Feb 2026 13:10:33 -0700 Subject: [PATCH 09/12] Coef pseudo time in namelist --- io/namelist_mod.F90 | 8 ++++++-- physics/fire_model_mod.F90 | 2 +- physics/level_set_mod.F90 | 7 +++---- tests/test7/namelist.fire | 1 + tests/test8/namelist.fire | 1 + tests/testx/namelist.fire | 1 + 6 files changed, 13 insertions(+), 7 deletions(-) diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index 8e71ef4..fc188ba 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -57,6 +57,7 @@ module namelist_mod integer :: fire_upwinding_reinit = 4 ! "numerical scheme (space) for reinitialization PDE: 1=WENO3, 2=WENO5, 3=hybrid WENO3-ENO1, 4=hybrid WENO5-ENO1" integer :: fire_lsm_band_ngp = 4 ! "number of grid points around lfn=0 that WENO5/3 is used (ENO1 elsewhere), ! for fire_upwinding_reinit=4,5 and fire_upwinding=8,9 options" + real :: reinit_pseudot_coef = 0.01 ! Coefficient for the pseudo time integer :: fast_dist_reinit_opt = 0 ! Fast distance reinitialization method (or eikonal solver): 0) None, 1) FSM integer :: fast_dist_reinit_freq = 600 ! Number of time steps to perform a reinit with fast distance reinit method @@ -184,6 +185,7 @@ subroutine Broadcast_nml (this) call Broadcast_integer (this%fire_lsm_reinit_iter) call Broadcast_integer (this%fire_upwinding_reinit) call Broadcast_integer (this%fire_lsm_band_ngp) + call Broadcast_real (this%reinit_pseudot_coef) call Broadcast_integer (this%fast_dist_reinit_opt) call Broadcast_integer (this%fast_dist_reinit_freq) call Broadcast_logical (this%fire_lsm_zcoupling) @@ -406,7 +408,7 @@ subroutine Init_fire_block (this, file_name) fast_dist_reinit_opt, fast_dist_reinit_freq, fire_viscosity_ngp, wind_vinterp_opt, hinterp_opt, ideal_opt, & fuel_opt, ros_opt, fmc_opt, emis_opt, fmoist_freq real :: fire_atm_feedback, fire_viscosity, fire_lsm_zcoupling_ref, fire_viscosity_bg, fire_viscosity_band, & - fmoist_dt, fire_wind_height, frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c + fmoist_dt, fire_wind_height, frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, reinit_pseudot_coef logical :: fire_lsm_reinit, fire_lsm_zcoupling, fmoist_run, fire_is_real_perim ! ignitions @@ -427,7 +429,7 @@ subroutine Init_fire_block (this, file_name) fire_lsm_band_ngp, fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_viscosity_bg, fire_viscosity_band, & fire_viscosity_ngp, fmoist_run, fmoist_freq, fmoist_dt, fire_wind_height, fire_is_real_perim, & frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, ideal_opt, fuel_opt, ros_opt, fmc_opt, emis_opt, & - wind_vinterp_opt, hinterp_opt, & + wind_vinterp_opt, hinterp_opt, reinit_pseudot_coef, & ! Ignitions fire_num_ignitions, & ! Ignition 1 @@ -459,6 +461,7 @@ subroutine Init_fire_block (this, file_name) fire_lsm_reinit_iter = this%fire_lsm_reinit_iter fire_upwinding_reinit = this%fire_upwinding_reinit fire_lsm_band_ngp = this%fire_lsm_band_ngp + reinit_pseudot_coef = this%reinit_pseudot_coef fast_dist_reinit_opt = this%fast_dist_reinit_opt fast_dist_reinit_freq = this%fast_dist_reinit_freq fire_lsm_zcoupling = this%fire_lsm_zcoupling @@ -553,6 +556,7 @@ subroutine Init_fire_block (this, file_name) this%fire_lsm_reinit_iter = fire_lsm_reinit_iter this%fire_upwinding_reinit = fire_upwinding_reinit this%fire_lsm_band_ngp = fire_lsm_band_ngp + this%reinit_pseudot_coef = reinit_pseudot_coef this%fast_dist_reinit_opt = fast_dist_reinit_opt this%fast_dist_reinit_freq = fast_dist_reinit_freq this%fire_lsm_zcoupling = fire_lsm_zcoupling diff --git a/physics/fire_model_mod.F90 b/physics/fire_model_mod.F90 index 4efa9d7..225bb84 100644 --- a/physics/fire_model_mod.F90 +++ b/physics/fire_model_mod.F90 @@ -113,7 +113,7 @@ subroutine Advance_fire_model (config_flags, grid) ifds, ifde, jfds, jfde, time_start, grid%dt, grid%dx, grid%dy, config_flags%fire_upwinding_reinit, & config_flags%fire_lsm_reinit_iter, config_flags%fire_lsm_band_ngp, grid%lfn, grid%lfn_2, grid%lfn_s0, & grid%lfn_s1, grid%lfn_s2, grid%lfn_s3, grid%lfn_out, grid%tign_g, grid%cart_comm, & - grid%ifps, grid%ifpe, grid%jfps, grid%jfpe) + grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, config_flags%reinit_pseudot_coef) if (DEBUG_LOCAL) call Print_message ('calling Copy_lfnout_to_lfn...') !$OMP PARALLEL DO & diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 94809e2..1e5bd67 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -618,7 +618,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if ifds, ifde, jfds, jfde, ts, dt, dx, dy, fire_upwinding_reinit, & fire_lsm_reinit_iter, fire_lsm_band_ngp, lfn_in, lfn_2, lfn_s0, & lfn_s1, lfn_s2, lfn_s3, lfn_out, tign, cart_comm, & - ifps, ifpe, jfps, jfpe) + ifps, ifpe, jfps, jfpe, reinit_pseudot_coef) ! Purpose: Level-set function reinitialization ! @@ -640,7 +640,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_2, lfn_s0, lfn_s1, lfn_s2, lfn_s3 real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_out - real, intent (in) :: dx, dy, ts, dt + real, intent (in) :: reinit_pseudot_coef, dx, dy, ts, dt real :: dt_s, threshold_hlu integer :: nts, i, j, ij, ifts, ifte, jfts, jfte @@ -683,8 +683,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if end do !$OMP END PARALLEL DO - dt_s = 0.01 * dx - dt_s = 0.0001 * dx + dt_s = reinit_pseudot_coef * dx ! iterate to solve to steady state reinit PDE ! 1 iter each time step is enoguh diff --git a/tests/test7/namelist.fire b/tests/test7/namelist.fire index d6931c7..7c02c76 100644 --- a/tests/test7/namelist.fire +++ b/tests/test7/namelist.fire @@ -37,5 +37,6 @@ fire_viscosity = 0.4, ! artificial viscosity in level set method (max 1, needed with fire_upwinding=0) fire_upwinding = 9, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian wind_vinterp_opt = 0, + reinit_pseudot_coef = 0.0001 / diff --git a/tests/test8/namelist.fire b/tests/test8/namelist.fire index 502787d..62d4b6f 100644 --- a/tests/test8/namelist.fire +++ b/tests/test8/namelist.fire @@ -39,5 +39,6 @@ fmoist_run = .true., ! Run fuel moisture model fmoist_freq = 1, ! Fuel moisture model updated every time step wind_vinterp_opt = 0 + reinit_pseudot_coef = 0.0001 / diff --git a/tests/testx/namelist.fire b/tests/testx/namelist.fire index 5dfd895..0346d66 100644 --- a/tests/testx/namelist.fire +++ b/tests/testx/namelist.fire @@ -37,6 +37,7 @@ fire_viscosity = 0.4, ! artificial viscosity in level set method (max 1, needed with fire_upwinding=0) fire_upwinding = 9, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian wind_vinterp_opt = 0 + reinit_pseudot_coef = 0.0001 / &ideal From cb5f8933feae2623e31e0efab4d6849e43ba8127 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Sun, 15 Feb 2026 15:11:05 -0700 Subject: [PATCH 10/12] Adds devel block to namelist --- io/namelist_mod.F90 | 45 ++++++++++++++++++++++++++++++++++++-- physics/fire_model_mod.F90 | 3 +-- 2 files changed, 44 insertions(+), 4 deletions(-) diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index fc188ba..622d27e 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -78,6 +78,7 @@ module namelist_mod real :: fmoist_dt = 600.0 ! moisture model time step [s] integer :: ideal_opt = 0 ! 0) real world, 1) ideal + integer :: devel_opt = 0 ! 0) Standard nml options, 1) reads options in the devel nml block ! Objects integer :: fuel_opt = FUEL_ANDERSON ! 1) Anderson 13 @@ -155,10 +156,14 @@ module namelist_mod ! Atm block integer :: interval_atm = -1 ! Time step [s] of the atm (or frequency to read atm data if offline) integer :: kde = 1 ! Number of atm vertical levels + + ! Devel block + integer :: check_isolated_neg_lfn = 0 ! 0) Nothing, 1) Check for isolated negative values of the level set function contains procedure, public :: Broadcast_nml => Broadcast_nml procedure, public :: Check_nml => Check_nml procedure, public :: Initialization => Init_namelist + procedure, public :: Init_devel_block => Init_devel_block procedure, public :: Init_fire_block => Init_fire_block procedure, public :: Init_ideal_block => Init_ideal_block procedure, public :: Init_time_block => Init_time_block @@ -204,6 +209,7 @@ subroutine Broadcast_nml (this) call Broadcast_real (this%fuelmc_c) call Broadcast_integer (this%ideal_opt) + call Broadcast_integer (this%devel_opt) call Broadcast_integer (this%fuel_opt) call Broadcast_integer (this%ros_opt) call Broadcast_integer (this%emis_opt) @@ -300,6 +306,8 @@ subroutine Broadcast_nml (this) call Broadcast_integer (this%interval_atm) call Broadcast_integer (this%kde) + ! Devel block + call Broadcast_integer (this%check_isolated_neg_lfn) contains subroutine Broadcast_integer (val) @@ -397,6 +405,36 @@ subroutine Init_atm_block (this, file_name) end subroutine Init_atm_block + subroutine Init_devel_block (this, file_name) + + implicit none + + class (namelist_t), intent (in out) :: this + character (len = *), intent (in) :: file_name + + integer :: check_isolated_neg_lfn + integer :: unit_nml, io_stat + character (len = :), allocatable :: msg + + namelist /devel/ check_isolated_neg_lfn + + + check_isolated_neg_lfn = this%check_isolated_neg_lfn + + open (newunit = unit_nml, file = trim (file_name), action = 'read', iostat = io_stat) + if (io_stat /= 0) then + msg = 'Problems opening namelist file ' // trim (file_name) + call Stop_simulation (msg) + end if + + read (unit_nml, nml = devel, iostat = io_stat) + if (io_stat /= 0) call Stop_simulation ('Problems reading namelist devel block') + close (unit_nml) + + this%check_isolated_neg_lfn = check_isolated_neg_lfn + + end subroutine Init_devel_block + subroutine Init_fire_block (this, file_name) implicit none @@ -405,7 +443,7 @@ subroutine Init_fire_block (this, file_name) character (len = *), intent (in) :: file_name integer :: fire_print_msg, fire_upwinding, fire_lsm_reinit_iter, fire_upwinding_reinit, fire_lsm_band_ngp, & - fast_dist_reinit_opt, fast_dist_reinit_freq, fire_viscosity_ngp, wind_vinterp_opt, hinterp_opt, ideal_opt, & + fast_dist_reinit_opt, fast_dist_reinit_freq, fire_viscosity_ngp, wind_vinterp_opt, hinterp_opt, ideal_opt, devel_opt, & fuel_opt, ros_opt, fmc_opt, emis_opt, fmoist_freq real :: fire_atm_feedback, fire_viscosity, fire_lsm_zcoupling_ref, fire_viscosity_bg, fire_viscosity_band, & fmoist_dt, fire_wind_height, frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, reinit_pseudot_coef @@ -428,7 +466,7 @@ subroutine Init_fire_block (this, file_name) fast_dist_reinit_opt, fast_dist_reinit_freq, fire_lsm_reinit_iter, fire_upwinding_reinit, & fire_lsm_band_ngp, fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_viscosity_bg, fire_viscosity_band, & fire_viscosity_ngp, fmoist_run, fmoist_freq, fmoist_dt, fire_wind_height, fire_is_real_perim, & - frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, ideal_opt, fuel_opt, ros_opt, fmc_opt, emis_opt, & + frac_fburnt_to_smoke, fuelmc_g, fuelmc_g_live, fuelmc_c, ideal_opt, devel_opt, fuel_opt, ros_opt, fmc_opt, emis_opt, & wind_vinterp_opt, hinterp_opt, reinit_pseudot_coef, & ! Ignitions fire_num_ignitions, & @@ -480,6 +518,7 @@ subroutine Init_fire_block (this, file_name) fuelmc_c = this%fuelmc_c ideal_opt = this%ideal_opt + devel_opt = this%devel_opt fuel_opt = this%fuel_opt ros_opt = this%ros_opt @@ -575,6 +614,7 @@ subroutine Init_fire_block (this, file_name) this%fuelmc_c = fuelmc_c this%ideal_opt = ideal_opt + this%devel_opt = devel_opt this%fuel_opt = fuel_opt this%ros_opt = ros_opt @@ -790,6 +830,7 @@ subroutine Init_namelist (this, file_name) call this%Init_fire_block (file_name = trim (file_name)) call this%Init_atm_block (file_name = trim (file_name)) if (this%ideal_opt > 0) call this%Init_ideal_block (file_name = trim (file_name)) + if (this%devel_opt > 0) call this%Init_devel_block (file_name = trim (file_name)) call this%Check_nml () diff --git a/physics/fire_model_mod.F90 b/physics/fire_model_mod.F90 index 225bb84..a44c8c2 100644 --- a/physics/fire_model_mod.F90 +++ b/physics/fire_model_mod.F90 @@ -29,7 +29,6 @@ subroutine Advance_fire_model (config_flags, grid) integer :: ij, ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme real :: tbound, time_start - integer, parameter :: CHECK_ISOLATED_NEG_LFN = 0 logical, parameter :: DEBUG_LOCAL = .false. @@ -132,7 +131,7 @@ subroutine Advance_fire_model (config_flags, grid) call Do_halo_exchange_with_corners (grid%lfn, ifms, ifme, jfms, jfme, grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, N_POINTS_IN_HALO, grid%cart_comm) #endif - if (CHECK_ISOLATED_NEG_LFN == 1) call Check_isolated_negative_lfn (grid) + if (config_flags%check_isolated_neg_lfn == 1) call Check_isolated_negative_lfn (grid) if (DEBUG_LOCAL) call Print_message ('calling Ignite_prescribed_fires...') !$OMP PARALLEL DO & From 14087568bb31b8c5fc1889e9a98ce1144ac84abc Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Sun, 15 Feb 2026 16:47:56 -0700 Subject: [PATCH 11/12] Enhanced output available --- io/namelist_mod.F90 | 8 ++++++-- physics/fire_model_mod.F90 | 4 ++-- physics/level_set_mod.F90 | 29 +++++++++++++++++------------ state/state_mod.F90 | 18 ++++++++++++++++++ 4 files changed, 43 insertions(+), 16 deletions(-) diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index 622d27e..9295e4a 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -159,6 +159,7 @@ module namelist_mod ! Devel block integer :: check_isolated_neg_lfn = 0 ! 0) Nothing, 1) Check for isolated negative values of the level set function + integer :: output_level = 0 ! 0) Standard output, >0) Specialized output contains procedure, public :: Broadcast_nml => Broadcast_nml procedure, public :: Check_nml => Check_nml @@ -308,6 +309,7 @@ subroutine Broadcast_nml (this) ! Devel block call Broadcast_integer (this%check_isolated_neg_lfn) + call Broadcast_integer (this%output_level) contains subroutine Broadcast_integer (val) @@ -412,14 +414,15 @@ subroutine Init_devel_block (this, file_name) class (namelist_t), intent (in out) :: this character (len = *), intent (in) :: file_name - integer :: check_isolated_neg_lfn + integer :: check_isolated_neg_lfn, output_level integer :: unit_nml, io_stat character (len = :), allocatable :: msg - namelist /devel/ check_isolated_neg_lfn + namelist /devel/ check_isolated_neg_lfn, output_level check_isolated_neg_lfn = this%check_isolated_neg_lfn + output_level = this%output_level open (newunit = unit_nml, file = trim (file_name), action = 'read', iostat = io_stat) if (io_stat /= 0) then @@ -432,6 +435,7 @@ subroutine Init_devel_block (this, file_name) close (unit_nml) this%check_isolated_neg_lfn = check_isolated_neg_lfn + this%output_level = output_level end subroutine Init_devel_block diff --git a/physics/fire_model_mod.F90 b/physics/fire_model_mod.F90 index a44c8c2..c6dfed2 100644 --- a/physics/fire_model_mod.F90 +++ b/physics/fire_model_mod.F90 @@ -52,7 +52,7 @@ subroutine Advance_fire_model (config_flags, grid) config_flags%fire_upwinding, config_flags%fire_viscosity, config_flags%fire_viscosity_bg, config_flags%fire_viscosity_band, & config_flags%fire_viscosity_ngp, config_flags%fire_lsm_band_ngp, tbound, grid%lfn, grid%lfn_0, grid%lfn_1, grid%lfn_2, & grid%lfn_out, grid%tign_g, grid%ros, grid%uf, grid%vf, grid%dzdxf, grid%dzdyf, grid%ros_param, grid%cart_comm, & - grid%ifps, grid%ifpe, grid%jfps, grid%jfpe) + grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, grid%grad_norm_ls) if (DEBUG_LOCAL) call Print_message ('calling Stop_if_close_to_bdy...') !$OMP PARALLEL DO & @@ -112,7 +112,7 @@ subroutine Advance_fire_model (config_flags, grid) ifds, ifde, jfds, jfde, time_start, grid%dt, grid%dx, grid%dy, config_flags%fire_upwinding_reinit, & config_flags%fire_lsm_reinit_iter, config_flags%fire_lsm_band_ngp, grid%lfn, grid%lfn_2, grid%lfn_s0, & grid%lfn_s1, grid%lfn_s2, grid%lfn_s3, grid%lfn_out, grid%tign_g, grid%cart_comm, & - grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, config_flags%reinit_pseudot_coef) + grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, config_flags%reinit_pseudot_coef, grid%grad_norm_reinit) if (DEBUG_LOCAL) call Print_message ('calling Copy_lfnout_to_lfn...') !$OMP PARALLEL DO & diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index 1e5bd67..d59e646 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -427,7 +427,7 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & num_tiles, i_start, i_end, j_start, j_end, ts, dt, dx, dy, fire_upwinding, fire_viscosity, & fire_viscosity_bg, fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, & tbound, lfn_in, lfn_0, lfn_1, lfn_2, lfn_out, tign, ros, uf, vf, dzdxf, dzdyf, ros_model, cart_comm, & - ifps, ifpe, jfps, jfpe) + ifps, ifpe, jfps, jfpe, grad_norm_ls) ! Purpose: Advance the level set function from time ts to time ts + dt @@ -439,7 +439,7 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & real, intent(in) :: fire_viscosity, fire_viscosity_bg, fire_viscosity_band real, dimension(ifms:ifme, jfms:jfme), intent (in) :: uf, vf, dzdxf, dzdyf real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign, lfn_1, lfn_2, lfn_0 - real, dimension(ifms:ifme, jfms:jfme), intent (out) :: lfn_out, ros + real, dimension(ifms:ifme, jfms:jfme), intent (out) :: lfn_out, ros, grad_norm_ls real, intent (in) :: dx, dy, ts, dt real, intent (out) :: tbound class (ros_t), intent (in) :: ros_model @@ -488,7 +488,7 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & ifms, ifme, jfms, jfme, ts, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) tbound_min = min(tbound_min, tbound_thread) end do @@ -531,7 +531,7 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & ifms,ifme,jfms,jfme, ts + dt, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) tbound_min = min(tbound_min, tbound_thread) end do @@ -573,7 +573,7 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds,ifde,jfds,jfde, ifts, ifte, jfts, jfte, & ifms, ifme, jfms, jfme, ts + dt, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) tbound_min = min(tbound_min, tbound_thread) end do @@ -618,7 +618,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if ifds, ifde, jfds, jfde, ts, dt, dx, dy, fire_upwinding_reinit, & fire_lsm_reinit_iter, fire_lsm_band_ngp, lfn_in, lfn_2, lfn_s0, & lfn_s1, lfn_s2, lfn_s3, lfn_out, tign, cart_comm, & - ifps, ifpe, jfps, jfpe, reinit_pseudot_coef) + ifps, ifpe, jfps, jfpe, reinit_pseudot_coef, grad_norm_reinit) ! Purpose: Level-set function reinitialization ! @@ -640,6 +640,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_2, lfn_s0, lfn_s1, lfn_s2, lfn_s3 real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_out + real, dimension (ifms:ifme, jfms:jfme), intent (out) :: grad_norm_reinit real, intent (in) :: reinit_pseudot_coef, dx, dy, ts, dt real :: dt_s, threshold_hlu @@ -700,7 +701,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & lfn_s0, lfn_s3, lfn_s3, lfn_s1, 1.0 / 3.0, & ! sign funcition, initial ls, current stage ls, next stage advanced ls, RK coefficient - fire_upwinding_reinit) + fire_upwinding_reinit, grad_norm_reinit) end do !$OMP END PARALLEL DO @@ -733,7 +734,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & lfn_s0, lfn_s3, lfn_s1, lfn_s2, 1.0 / 2.0, & - fire_upwinding_reinit) + fire_upwinding_reinit, grad_norm_reinit) end do !$OMP END PARALLEL DO @@ -766,7 +767,7 @@ subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, if call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & lfn_s0, lfn_s3, lfn_s2, lfn_s3, 1.0, & - fire_upwinding_reinit) + fire_upwinding_reinit, grad_norm_reinit) end do !$OMP END PARALLEL DO @@ -811,7 +812,7 @@ end subroutine Reinit_level_set subroutine Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, lfn_s0, & - lfn_ini, lfn_curr, lfn_fin, rk_coeff, fire_upwinding_reinit) + lfn_ini, lfn_curr, lfn_fin, rk_coeff, fire_upwinding_reinit, grad_norm_reinit) ! Calculates right-hand-side forcing and advances a RK-stage the level-set reinitialization PDE @@ -822,6 +823,7 @@ subroutine Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & integer, intent (in) :: fire_upwinding_reinit real, dimension (ifms:ifme, jfms:jfme), intent (in) :: lfn_s0, lfn_ini, lfn_curr real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_fin + real, dimension (ifms:ifme, jfms:jfme), intent (out) :: grad_norm_reinit real, intent (in) :: dx, dy, dt_s, threshold_hlu, rk_coeff integer :: i, j @@ -930,6 +932,7 @@ subroutine Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & end select end if grad = sqrt (diff2x * diff2x + diff2y * diff2y) + grad_norm_reinit(i, j) = grad tend_r = lfn_s0(i, j) * (1.0 - grad) lfn_fin(i, j) = lfn_ini(i, j) + (dt_s * rk_coeff) * tend_r end do @@ -969,7 +972,7 @@ end subroutine Update_ignition_times subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfms, jfme, & t, dt, dx, dy, fire_upwinding, fire_viscosity, fire_viscosity_bg, & - fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, lfn, tbound, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, lfn, tbound, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) ! compute the right hand side of the level set equation @@ -980,7 +983,7 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm real, intent (in) :: fire_viscosity, fire_viscosity_bg, fire_viscosity_band, t, dt, dx, dy real, dimension(ifms:ifme, jfms:jfme), intent (in) :: uf, vf, dzdxf, dzdyf real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: lfn - real, dimension(ifms:ifme, jfms:jfme), intent (out) :: tend, ros + real, dimension(ifms:ifme, jfms:jfme), intent (out) :: tend, ros, grad_norm_ls real, intent (out) :: tbound class (ros_t), intent (in) :: ros_model @@ -1177,6 +1180,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm end select end if + grad_norm_ls(i, j) = grad + ! Calc normal scale = sqrt (grad ** 2.0 + EPS) nvx = diff2x / scale diff --git a/state/state_mod.F90 b/state/state_mod.F90 index c0f5e61..eb9811b 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -73,6 +73,8 @@ module state_mod real, dimension(:, :), allocatable :: nfuel_cat ! "fuel data" real, dimension(:, :), allocatable :: fuel_time ! "fuel" real, dimension(:, :), allocatable :: emis_smoke + real, dimension(:, :), allocatable :: grad_norm_ls ! Gracient norm of the level set function used to propagate level set function + real, dimension(:, :), allocatable :: grad_norm_reinit ! Gracient norm of the level set function used to reinitilize the level set function class (fuel_t), allocatable :: fuels class (ros_t), allocatable :: ros_param @@ -98,6 +100,9 @@ module state_mod integer :: ny ! "number of latitudinal grid points" "1" real :: cen_lat, cen_lon + ! Output + integer :: output_level + ! For MPI tasks integer :: cart_comm ! The MPI communicator with the domain decomposition integer :: ntasks ! Number of MPI tasks @@ -173,6 +178,8 @@ subroutine Allocate_vars (this, ifms, ifme, jfms, jfme) allocate (this%dzdyf(ifms:ifme, jfms:jfme)) allocate (this%nfuel_cat(ifms:ifme, jfms:jfme)) allocate (this%emis_smoke(ifms:ifme, jfms:jfme)) + allocate (this%grad_norm_ls(ifms:ifme, jfms:jfme)) + allocate (this%grad_norm_reinit(ifms:ifme, jfms:jfme)) end subroutine Allocate_vars @@ -594,6 +601,9 @@ subroutine Init_domain (this, config_flags, geogrid, & if (DEBUG_LOCAL) call Print_message (' Setting clock...') call this%Set_time_stamps (config_flags) + ! Output + this%output_level = config_flags%output_level + if (DEBUG_LOCAL) call this%Print() if (DEBUG_LOCAL) call Print_message ('Leaving Init_domain...') @@ -944,6 +954,14 @@ subroutine Save_state (this) call Add_netcdf_var_mpi (file_output, this%nx, this%ny, this%ifps, this%ifpe, this%jfps, this%jfpe, 'nfuel_cat', & this%nfuel_cat(this%ifps:this%ifpe, this%jfps:this%jfpe)) + if (this%output_level > 0) then + call Add_netcdf_var_mpi (file_output, this%nx, this%ny, this%ifps, this%ifpe, this%jfps, this%jfpe, 'grad_norm_ls', & + this%grad_norm_ls(this%ifps:this%ifpe, this%jfps:this%jfpe)) + + call Add_netcdf_var_mpi (file_output, this%nx, this%ny, this%ifps, this%ifpe, this%jfps, this%jfpe, 'grad_norm_reinit', & + this%grad_norm_reinit(this%ifps:this%ifpe, this%jfps:this%jfpe)) + end if + if (DEBUG_LOCAL) call Print_message ('Leaving Save_state...') end subroutine Save_state From 757150f24d3448355921b4b343bc056dc3a77bb4 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Wed, 18 Feb 2026 14:01:24 -0700 Subject: [PATCH 12/12] Calc bulk errors in the level set function --- physics/fire_model_mod.F90 | 3 +- physics/level_set_mod.F90 | 161 ++++++++++++++++++++++++++++++++----- share/mpi_mod.F90 | 25 +++++- state/state_mod.F90 | 5 ++ 4 files changed, 174 insertions(+), 20 deletions(-) diff --git a/physics/fire_model_mod.F90 b/physics/fire_model_mod.F90 index c6dfed2..0e5679f 100644 --- a/physics/fire_model_mod.F90 +++ b/physics/fire_model_mod.F90 @@ -52,7 +52,8 @@ subroutine Advance_fire_model (config_flags, grid) config_flags%fire_upwinding, config_flags%fire_viscosity, config_flags%fire_viscosity_bg, config_flags%fire_viscosity_band, & config_flags%fire_viscosity_ngp, config_flags%fire_lsm_band_ngp, tbound, grid%lfn, grid%lfn_0, grid%lfn_1, grid%lfn_2, & grid%lfn_out, grid%tign_g, grid%ros, grid%uf, grid%vf, grid%dzdxf, grid%dzdyf, grid%ros_param, grid%cart_comm, & - grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, grid%grad_norm_ls) + grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, grid%grad_norm_ls, grid%grad_norm_residual_sq_sum, & + grid%grad_norm_residual_sq_sum_band, grid%grad_norm_residual_rms_band) if (DEBUG_LOCAL) call Print_message ('calling Stop_if_close_to_bdy...') !$OMP PARALLEL DO & diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index d59e646..5aa7bdf 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -21,7 +21,7 @@ module level_set_mod use constants_mod, only : PI #ifdef DM_PARALLEL - use mpi_mod, only : Do_halo_exchange, Sum_across_mpi_tasks, Max_across_mpi_tasks + use mpi_mod, only : Do_halo_exchange, Sum_across_mpi_tasks, Max_across_mpi_tasks, Min_across_mpi_tasks #endif implicit none @@ -427,7 +427,8 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & num_tiles, i_start, i_end, j_start, j_end, ts, dt, dx, dy, fire_upwinding, fire_viscosity, & fire_viscosity_bg, fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, & tbound, lfn_in, lfn_0, lfn_1, lfn_2, lfn_out, tign, ros, uf, vf, dzdxf, dzdyf, ros_model, cart_comm, & - ifps, ifpe, jfps, jfpe, grad_norm_ls) + ifps, ifpe, jfps, jfpe, grad_norm_ls, grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, & + grad_norm_residual_rms_band) ! Purpose: Advance the level set function from time ts to time ts + dt @@ -441,15 +442,16 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign, lfn_1, lfn_2, lfn_0 real, dimension(ifms:ifme, jfms:jfme), intent (out) :: lfn_out, ros, grad_norm_ls real, intent (in) :: dx, dy, ts, dt - real, intent (out) :: tbound + real, intent (out) :: tbound, grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band class (ros_t), intent (in) :: ros_model ! to store tendency (rhs of the level set pde) real, dimension(ifms:ifme, jfms:jfme) :: tend - real :: tbound2, tbound3, tbound_thread, tbound_min - integer :: i, j, ij, ifts, ifte, jfts, jfte + real :: tbound2, tbound3, tbound_thread, tbound_min, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, & + np_band, gmin, gsum + integer :: i, j, ij, ifts, ifte, jfts, jfte, np_band_local character (len = 128) :: msg - logical, parameter :: DEBUG_LOCAL = .false. + logical, parameter :: DEBUG_LOCAL = .false., PRINT_ERRORS = .false. if (DEBUG_LOCAL) call Print_message ('Entering sub Prop_level_set...') @@ -476,9 +478,15 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & ! Runge-Kutta step 1 if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 1...') - tbound_min = huge(tbound_min) + tbound_min = huge (tbound_min) + grad_norm_residual_sq_sum = 0.0 + grad_norm_residual_sq_sum_band = 0.0 + np_band = 0.0 !$OMP PARALLEL DO & - !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread, grad_norm_residual_sq_sum_local, & + !$OMP grad_norm_residual_sq_sum_band_local, np_band_local) & + !$OMP REDUCTION(MIN: tbound_min) & + !$OMP REDUCTION(+: grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, np_band) do ij = 1, num_tiles ifts = i_start(ij) ifte = i_end(ij) @@ -488,12 +496,44 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & ifms, ifme, jfms, jfme, ts, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) tbound_min = min(tbound_min, tbound_thread) + grad_norm_residual_sq_sum = grad_norm_residual_sq_sum + grad_norm_residual_sq_sum_local + grad_norm_residual_sq_sum_band = grad_norm_residual_sq_sum_band + grad_norm_residual_sq_sum_band_local + np_band = np_band + real (np_band_local) end do !$OMP END PARALLEL DO + +#ifdef DM_PARALLEL + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum, cart_comm, gsum) + grad_norm_residual_sq_sum = gsum + + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum_band, cart_comm, gsum) + grad_norm_residual_sq_sum_band = gsum + + call Sum_across_mpi_tasks (np_band, cart_comm, gsum) + np_band = gsum + + call Min_across_mpi_tasks (tbound_min, cart_comm, gmin) + tbound_min = gmin +#endif + tbound = tbound_min + grad_norm_residual_sq_sum = sqrt (grad_norm_residual_sq_sum * dx * dy) + grad_norm_residual_sq_sum_band = sqrt (grad_norm_residual_sq_sum_band * dx * dy) + if (np_band > 0.0) then + grad_norm_residual_rms_band = sqrt (grad_norm_residual_sq_sum_band / np_band) + else + grad_norm_residual_rms_band = 0.0 + end if + + if (PRINT_ERRORS) then + write (msg, '(a, 3f12.6)') 'grad_norm_residual_sq_sum rk1, band, rms_band = ', grad_norm_residual_sq_sum, & + grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band + call Print_message (msg) + end if !$OMP PARALLEL DO & !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) @@ -519,9 +559,15 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & ! Runge-Kutta step 2 if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 2...') - tbound_min = huge(tbound_min) + tbound_min = huge (tbound_min) + grad_norm_residual_sq_sum = 0.0 + grad_norm_residual_sq_sum_band = 0.0 + np_band = 0.0 !$OMP PARALLEL DO & - !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread, grad_norm_residual_sq_sum_local, & + !$OMP grad_norm_residual_sq_sum_band_local, np_band_local) & + !$OMP REDUCTION(MIN: tbound_min) & + !$OMP REDUCTION(+: grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, np_band) do ij = 1, num_tiles ifts = i_start(ij) ifte = i_end(ij) @@ -531,12 +577,44 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & ifms,ifme,jfms,jfme, ts + dt, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) - tbound_min = min(tbound_min, tbound_thread) + tbound_min = min(tbound_min, tbound_thread) + grad_norm_residual_sq_sum = grad_norm_residual_sq_sum + grad_norm_residual_sq_sum_local + grad_norm_residual_sq_sum_band = grad_norm_residual_sq_sum_band + grad_norm_residual_sq_sum_band_local + np_band = np_band + real (np_band_local) end do !$OMP END PARALLEL DO + +#ifdef DM_PARALLEL + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum, cart_comm, gsum) + grad_norm_residual_sq_sum = gsum + + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum_band, cart_comm, gsum) + grad_norm_residual_sq_sum_band = gsum + + call Sum_across_mpi_tasks (np_band, cart_comm, gsum) + np_band = gsum + + call Min_across_mpi_tasks (tbound_min, cart_comm, gmin) + tbound_min = gmin +#endif + tbound2 = tbound_min + grad_norm_residual_sq_sum = sqrt (grad_norm_residual_sq_sum * dx * dy) + grad_norm_residual_sq_sum_band = sqrt (grad_norm_residual_sq_sum_band * dx * dy) + if (np_band > 0.0) then + grad_norm_residual_rms_band = sqrt (grad_norm_residual_sq_sum_band / np_band) + else + grad_norm_residual_rms_band = 0.0 + end if + + if (PRINT_ERRORS) then + write (msg, '(a, 3f12.6)') 'grad_norm_residual_sq_sum rk2, band, rms_band = ', grad_norm_residual_sq_sum, & + grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band + call Print_message (msg) + end if !$OMP PARALLEL DO & !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) @@ -561,9 +639,15 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & ! Runge-Kutta step 3 if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 3...') - tbound_min = huge(tbound_min) + tbound_min = huge (tbound_min) + grad_norm_residual_sq_sum = 0.0 + grad_norm_residual_sq_sum_band = 0.0 + np_band = 0.0 !$OMP PARALLEL DO & - !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread, grad_norm_residual_sq_sum_local, & + !$OMP grad_norm_residual_sq_sum_band_local, np_band_local) & + !$OMP REDUCTION(MIN: tbound_min) & + !$OMP REDUCTION(+: grad_norm_residual_sq_sum, grad_norm_residual_sq_sum_band, np_band) do ij = 1, num_tiles ifts = i_start(ij) ifte = i_end(ij) @@ -573,12 +657,43 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & call Calc_tend_ls (ifds,ifde,jfds,jfde, ifts, ifte, jfts, jfte, & ifms, ifme, jfms, jfme, ts + dt, dt, dx, dy, fire_upwinding, & fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) tbound_min = min(tbound_min, tbound_thread) + grad_norm_residual_sq_sum = grad_norm_residual_sq_sum + grad_norm_residual_sq_sum_local + grad_norm_residual_sq_sum_band = grad_norm_residual_sq_sum_band + grad_norm_residual_sq_sum_band_local + np_band = np_band + real (np_band_local) end do !$OMP END PARALLEL DO + +#ifdef DM_PARALLEL + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum, cart_comm, gsum) + grad_norm_residual_sq_sum = gsum + + call Sum_across_mpi_tasks (grad_norm_residual_sq_sum_band, cart_comm, gsum) + grad_norm_residual_sq_sum_band = gsum + + call Sum_across_mpi_tasks (np_band, cart_comm, gsum) + np_band = gsum + + call Min_across_mpi_tasks (tbound_min, cart_comm, gmin) + tbound_min = gmin +#endif + tbound3 = tbound_min + grad_norm_residual_sq_sum = sqrt (grad_norm_residual_sq_sum * dx * dy) + grad_norm_residual_sq_sum_band = sqrt (grad_norm_residual_sq_sum_band * dx * dy) + if (np_band > 0.0) then + grad_norm_residual_rms_band = sqrt (grad_norm_residual_sq_sum_band / np_band) + else + grad_norm_residual_rms_band = 0.0 + end if + if (PRINT_ERRORS) then + write (msg, '(a, 3f12.6)') 'grad_norm_residual_sq_sum rk3, band, rms_band = ', grad_norm_residual_sq_sum, & + grad_norm_residual_sq_sum_band, grad_norm_residual_rms_band + call Print_message (msg) + end if !$OMP PARALLEL DO & !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) @@ -972,7 +1087,8 @@ end subroutine Update_ignition_times subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfms, jfme, & t, dt, dx, dy, fire_upwinding, fire_viscosity, fire_viscosity_bg, & - fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, lfn, tbound, tend, ros, uf, vf, dzdxf, dzdyf, ros_model, grad_norm_ls) + fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, lfn, tbound, tend, ros, uf, vf, dzdxf, dzdyf, & + ros_model, grad_norm_ls, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local, np_band_local) ! compute the right hand side of the level set equation @@ -980,11 +1096,12 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm integer, intent (in) :: ifms, ifme, jfms, jfme, its, ite, jts, jte, ids, ide, jds, jde, & fire_upwinding,fire_viscosity_ngp, fire_lsm_band_ngp + integer, intent (out) :: np_band_local real, intent (in) :: fire_viscosity, fire_viscosity_bg, fire_viscosity_band, t, dt, dx, dy real, dimension(ifms:ifme, jfms:jfme), intent (in) :: uf, vf, dzdxf, dzdyf real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: lfn real, dimension(ifms:ifme, jfms:jfme), intent (out) :: tend, ros, grad_norm_ls - real, intent (out) :: tbound + real, intent (out) :: tbound, grad_norm_residual_sq_sum_local, grad_norm_residual_sq_sum_band_local class (ros_t), intent (in) :: ros_model real, parameter :: EPS = epsilon (0.0), TOL = 100.0 * EPS @@ -1010,6 +1127,9 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm if (DEBUG_LOCAL) call Print_message ('starting ij loops') tbound = 0.0 + grad_norm_residual_sq_sum_local = 0.0 + grad_norm_residual_sq_sum_band_local = 0.0 + np_band_local = 0 do j = jts, jte do i = its, ite ! one sided differences @@ -1181,6 +1301,11 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm end if grad_norm_ls(i, j) = grad + grad_norm_residual_sq_sum_local = grad_norm_residual_sq_sum_local + (grad - 1.0) ** 2 + if (abs(lfn(i, j)) < threshold_hlu) then + grad_norm_residual_sq_sum_band_local = grad_norm_residual_sq_sum_band_local + (grad - 1.0) ** 2 + np_band_local = np_band_local + 1 + end if ! Calc normal scale = sqrt (grad ** 2.0 + EPS) diff --git a/share/mpi_mod.F90 b/share/mpi_mod.F90 index 27aa26e..b8a1280 100644 --- a/share/mpi_mod.F90 +++ b/share/mpi_mod.F90 @@ -7,7 +7,7 @@ module mpi_mod private public :: Calc_tasks_in_x_and_y, Calc_patch_dims, Gather_var2d, Do_halo_exchange, Do_halo_exchange_with_corners, & - Max_across_mpi_tasks, Sum_across_mpi_tasks, Distribute_var2d + Max_across_mpi_tasks, Sum_across_mpi_tasks, Distribute_var2d, Min_across_mpi_tasks contains @@ -538,6 +538,29 @@ subroutine Max_across_mpi_tasks (local_max, cart_comm, global_max) end subroutine Max_across_mpi_tasks + subroutine Min_across_mpi_tasks (local_min, cart_comm, global_min) + +#ifdef DM_PARALLEL + use mpi +#endif + + implicit none + + real, intent (in) :: local_min + integer, intent (in) :: cart_comm + real, intent (out) :: global_min + + integer :: ierr + + +#ifdef DM_PARALLEL + call MPI_Allreduce (local_min, global_min, 1, MPI_REAL, MPI_MIN, cart_comm, ierr) +#else + global_min = local_min +#endif + + end subroutine Min_across_mpi_tasks + subroutine Sum_across_mpi_tasks (local_sum, cart_comm, global_sum) #ifdef DM_PARALLEL diff --git a/state/state_mod.F90 b/state/state_mod.F90 index eb9811b..6f574f8 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -100,6 +100,11 @@ module state_mod integer :: ny ! "number of latitudinal grid points" "1" real :: cen_lat, cen_lon + ! Performance stats + real :: grad_norm_residual_sq_sum + real :: grad_norm_residual_sq_sum_band + real :: grad_norm_residual_rms_band + ! Output integer :: output_level