diff --git a/src/ebfm/core/LOOP_SNOW.py b/src/ebfm/core/LOOP_SNOW.py index 016c75c8..21c07416 100644 --- a/src/ebfm/core/LOOP_SNOW.py +++ b/src/ebfm/core/LOOP_SNOW.py @@ -7,6 +7,7 @@ from .compute_backend import get_backend, ComputeBackend from .constants import SECONDS_PER_HOUR +from .LOOP_SNOW_common import _DT_SAFETY_FACTOR from .LOOP_SNOW_kernels import ( _compaction_kernel, _heat_conduction_kernel, @@ -546,7 +547,7 @@ def heat_conduction(gpu=None): # Stability time step (CFL condition) # Layer 0: surface ghost layer, excluded from CFL condition dt_stab: NDArray[np.float64] = ( - 0.5 + _DT_SAFETY_FACTOR * np.min(c_eff[:, 1:], axis=1) * np.min(OUT["subZ"][:, 1:], axis=1) ** 2 / np.max(kk[:, 1:], axis=1) diff --git a/src/ebfm/core/LOOP_SNOW_common.py b/src/ebfm/core/LOOP_SNOW_common.py new file mode 100644 index 00000000..7b92053a --- /dev/null +++ b/src/ebfm/core/LOOP_SNOW_common.py @@ -0,0 +1,5 @@ +# SPDX-FileCopyrightText: 2026 EBFM Authors +# +# SPDX-License-Identifier: BSD-3-Clause + +_DT_SAFETY_FACTOR = 0.5 diff --git a/src/ebfm/core/LOOP_SNOW_gpu_kernels.py b/src/ebfm/core/LOOP_SNOW_gpu_kernels.py index 0117c405..446f4260 100644 --- a/src/ebfm/core/LOOP_SNOW_gpu_kernels.py +++ b/src/ebfm/core/LOOP_SNOW_gpu_kernels.py @@ -41,6 +41,8 @@ from .compute_backend import cuda +from .LOOP_SNOW_common import _DT_SAFETY_FACTOR + # Threads per block used for all 1-D column-parallel launches. _TPB = 128 @@ -509,7 +511,7 @@ def _heat_conduction_prep_kernel_gpu( # denom_bottom = c_eff[-1] * (0.25*subZ[-2] + 0.75*subZ[-1]) denom_bottom[i] = c_eff[i, nl - 1] * (0.25 * subZ[i, nl - 2] + 0.75 * subZ[i, nl - 1]) - # dt_stab = 0.5 * min(c_eff[1:]) * min(subZ[1:])**2 / max(kk[1:]) / dayseconds + # dt_stab = _DT_SAFETY_FACTOR * min(c_eff[1:]) * min(subZ[1:])**2 / max(kk[1:]) / dayseconds min_ceff = math.inf min_sz = math.inf max_kk = 0.0 @@ -520,7 +522,7 @@ def _heat_conduction_prep_kernel_gpu( min_sz = subZ[i, k] if kk[i, k] > max_kk: max_kk = kk[i, k] - dt_stab[i] = 0.5 * min_ceff * (min_sz * min_sz) / max_kk / dayseconds + dt_stab[i] = _DT_SAFETY_FACTOR * min_ceff * (min_sz * min_sz) / max_kk / dayseconds @cuda.jit diff --git a/src/ebfm/core/LOOP_SNOW_kernels.py b/src/ebfm/core/LOOP_SNOW_kernels.py index 4f792473..cb874a3a 100644 --- a/src/ebfm/core/LOOP_SNOW_kernels.py +++ b/src/ebfm/core/LOOP_SNOW_kernels.py @@ -22,6 +22,8 @@ from .compute_backend import njit, prange +from .LOOP_SNOW_common import _DT_SAFETY_FACTOR + @njit(parallel=True, cache=True) def _compaction_kernel( @@ -224,7 +226,7 @@ def _heat_conduction_prep_kernel( # denom_bottom = c_eff[-1] * (0.25*subZ[-2] + 0.75*subZ[-1]) denom_bottom[i] = c_eff[i, nl - 1] * (0.25 * subZ[i, nl - 2] + 0.75 * subZ[i, nl - 1]) - # dt_stab = 0.5 * min(c_eff[1:]) * min(subZ[1:])**2 / max(kk[1:]) / dayseconds + # dt_stab = _DT_SAFETY_FACTOR * min(c_eff[1:]) * min(subZ[1:])**2 / max(kk[1:]) / dayseconds min_ceff = math.inf min_sz = math.inf max_kk = 0.0 @@ -235,7 +237,7 @@ def _heat_conduction_prep_kernel( min_sz = subZ[i, k] if kk[i, k] > max_kk: max_kk = kk[i, k] - dt_stab[i] = 0.5 * min_ceff * (min_sz * min_sz) / max_kk / dayseconds + dt_stab[i] = _DT_SAFETY_FACTOR * min_ceff * (min_sz * min_sz) / max_kk / dayseconds @njit(parallel=True, cache=True)