From 1ca98b6525a79684d5f0121af95bead6f70b78a9 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Fri, 23 Oct 2026 14:41:11 -0400 Subject: [PATCH 1/4] Add Arrhenius temperature-dependent shear viscosity (fluid_pp%visc_model) fluid_pp(i)%visc_model = 1 gives fluid i a shear viscosity mu_i(T) = exp(C_i + D_i/T) (visc_c, visc_d), evaluated at the thermal-equilibrium mixture temperature from f_mixture_temperature, the closure heat conduction uses. It overrides Re_L(1)/Re_R(1) on the HLL and HLLC (5- and 6-equation) Riemann states and Re(1) in the time-step and stability-criterion paths. The exponent is capped at 50. The validator requires viscous, model_eqns 2/3, riemann_solver 1/2, cv > 0 and nonzero visc_d, and rejects non_newtonian, cyl_coord, bubbles_euler, ib, igr and chemistry. Adds golden test 6D29F444 (Viscous -> visc_model=1) and validator unit tests. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VtsGVvCfsDJz5dmtSk9dKt --- docs/documentation/case.md | 5 + docs/module_categories.json | 1 + src/common/m_constants.fpp | 2 + src/common/m_derived_types.fpp | 3 + src/common/m_global_parameters_common.fpp | 3 +- src/post_process/m_global_parameters.fpp | 3 + src/pre_process/m_global_parameters.fpp | 3 + src/simulation/m_arrhenius_viscosity.fpp | 52 +++++ src/simulation/m_data_output.fpp | 5 + src/simulation/m_global_parameters.fpp | 24 +++ src/simulation/m_riemann_solver_hll.fpp | 5 + src/simulation/m_riemann_solver_hllc.fpp | 9 + src/simulation/m_time_steppers.fpp | 3 + tests/6D29F444/golden-metadata.txt | 196 ++++++++++++++++++ tests/6D29F444/golden.txt | 16 ++ toolchain/mfc/case_validator.py | 68 ++++++ toolchain/mfc/params/definitions.py | 9 + toolchain/mfc/params/descriptions.py | 3 + .../mfc/params_tests/test_eos_families.py | 3 + toolchain/mfc/test/cases.py | 4 + toolchain/mfc/test_case_validator.py | 37 ++++ 21 files changed, 453 insertions(+), 1 deletion(-) create mode 100644 src/simulation/m_arrhenius_viscosity.fpp create mode 100644 tests/6D29F444/golden-metadata.txt create mode 100644 tests/6D29F444/golden.txt diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 259b0bed47..5135c7aa89 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -489,6 +489,9 @@ A particle cloud is a compact specification of a bed of identical circular (2D) | `Re(1)` * | Real | Shear viscosity of fluid. | | `Re(2)` * | Real | Volume viscosity of fluid. | | `k_therm` | Real | Thermal conductivity of fluid (Fourier heat conduction). | +| `visc_model` | Integer | Shear viscosity model of fluid (0 constant, 1 Arrhenius \f$\mu(T)\f$). | +| `visc_c` | Real | Arrhenius viscosity coefficient \f$C\f$ of fluid. | +| `visc_d` | Real | Arrhenius viscosity coefficient \f$D\f$ of fluid. | | `cv` ** | Real | Sffened-gas parameter $c_v$ of fluid. | | `qv` ** | Real | Stiffened-gas parameter $q$ of fluid. | | `qvp` ** | Real | Stiffened-gas parameter $q'$ of fluid. | @@ -513,6 +516,8 @@ Details of implementation of viscosity in MFC can be found in \cite Coralic15. - `fluid_pp(i)%%k_therm` sets the thermal conductivity of the $i$-th fluid, in units consistent with the rest of the (non-dimensional) case. A positive value on any fluid activates Fourier heat conduction, which adds \f$\nabla\cdot(k\nabla T)\f$ to the energy equation using the thermal-equilibrium mixture temperature and \f$k = \sum_i \alpha_i k_i\f$ (see @ref equations "Equations"). It requires `fluid_pp(i)%%cv` to be positive on every fluid that sets it (the mixture temperature is undefined without \f$c_v\f$), `model_eqns = 2` or `model_eqns = 3` (the mixture conductivity is weighted by volume fractions that `model_eqns = 1` does not carry), and `fluid_pp(i)%%eos` to be the stiffened-gas or ideal-gas equation of state. Heat conduction is independent of `viscous`: it can be enabled in an otherwise inviscid run. It is not supported with `igr`, nor with `chemistry` (which already carries its own mixture-averaged conduction through `chem_params%%diffusion`). +- `fluid_pp(i)%%visc_model = 1` makes the shear viscosity of the $i$-th fluid temperature dependent, \f$\mu_i(T) = \exp(C_i + D_i/T)\f$ with \f$C_i\f$ = `fluid_pp(i)%%visc_c` and \f$D_i\f$ = `fluid_pp(i)%%visc_d` (nonzero), evaluated at the thermal-equilibrium mixture temperature that heat conduction uses. It replaces \f$1/\mathrm{Re}_i(1)\f$ in the mixture \f$1/\mathrm{Re} = \sum_i \alpha_i \mu_i\f$ on the left and right Riemann states, so the face viscosity follows the reconstructed temperature, and the viscous time-step limit sees the cell value. `fluid_pp(i)%%Re(1)` must still be set (`viscous` requires it, as for `non_newtonian`) but its value is unused for this fluid; `Re(2)` (bulk viscosity) stays constant. It requires `viscous = T`, `model_eqns = 2` or `3`, `riemann_solver = 1` or `2`, `fluid_pp(i)%%cv > 0` on every fluid, and the stiffened-gas or ideal-gas equation of state. It is not supported with `non_newtonian`, `cyl_coord`, `bubbles_euler`, `ib`, `igr`, or `chemistry`. The exponent \f$C_i + D_i/T\f$ is capped at 50 so that \f$\mu\f$ stays finite as \f$T \to 0\f$. At an isothermal wall (`isothermal_in`/`isothermal_out`) the wall-face viscosity is evaluated from the reflected ghost state, i.e. the adjacent interior temperature, not from `Twall`. + - `fluid_pp(i)%%cv`, `fluid_pp(i)%%qv`, and `fluid_pp(i)%%qvp` define $c_v$, $q$, and $q'$ as parameters of $i$-th fluid that are used in stiffened gas equation of state. - `fluid_pp(i)%%G` is required for `hypoelasticity`. diff --git a/docs/module_categories.json b/docs/module_categories.json index 0a7dafb33b..14040ed3f0 100644 --- a/docs/module_categories.json +++ b/docs/module_categories.json @@ -23,6 +23,7 @@ "m_eos", "m_viscous", "m_hb_function", + "m_arrhenius_viscosity", "m_surface_tension", "m_reactive_burn", "m_bubbles", diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 83507ac387..f0b677b9c2 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -52,6 +52,8 @@ module m_constants !> Spatial support width of acoustic source, used in s_source_spatial real(wp), parameter :: acoustic_spatial_support_width = 2.5_wp real(wp), parameter :: dflt_vcfl_dt = 100._wp !< value of vcfl_dt when viscosity is off for computing adaptive timestep size + !> Cap on the Arrhenius viscosity exponent C + D/T, so exp() stays finite in single precision (log(huge) ~ 88) as T -> 0 + real(wp), parameter :: visc_exp_arg_cap = 50._wp !> The constant to scale the spectral level at the lower frequency bound real(wp), parameter :: broadband_spectral_level_constant = 20._wp !> The spectral level constant to correct the magnitude at each frequency to ensure the source is overall broadband diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 1232e3f5ae..f29f72db67 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -447,6 +447,9 @@ module m_derived_types real(wp) :: mu_min !< Lower viscosity clamp (inactive sentinel = dflt_real) real(wp) :: mu_max !< Upper viscosity clamp (required when non_newtonian) real(wp) :: mu_bulk !< Bulk viscosity for NN (inactive sentinel = dflt_real) + integer :: visc_model !< Shear viscosity model selector (visc_model_* in m_constants) + real(wp) :: visc_c !< Arrhenius mu(T) = exp(visc_c + visc_d/T): coefficient C + real(wp) :: visc_d !< Arrhenius mu(T) = exp(visc_c + visc_d/T): coefficient D end type physical_parameters !> Derived type annexing the physical parameters required for sub-grid bubble models diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 80e3c2923a..a0727a5a8c 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -16,7 +16,8 @@ module m_global_parameters_common use m_derived_types use m_thermochem, only: num_species use m_constants, only: num_fluids_max, model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, recon_type_weno, & - & recon_type_muscl, name_len, dflt_int, dflt_real, eos_stiffened_gas, eos_ideal_gas + & recon_type_muscl, name_len, dflt_int, dflt_real, eos_stiffened_gas, eos_ideal_gas, visc_model_constant, & + & visc_model_arrhenius implicit none diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index bd90c28ed8..e35b97e9f9 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -245,6 +245,9 @@ contains fluid_pp(i)%mu_min = dflt_real fluid_pp(i)%mu_max = dflt_real fluid_pp(i)%mu_bulk = dflt_real + fluid_pp(i)%visc_model = visc_model_constant + fluid_pp(i)%visc_c = dflt_real + fluid_pp(i)%visc_d = dflt_real end do num_particle_clouds = 0 diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index b358516fbf..5b7396b9ef 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -462,6 +462,9 @@ contains fluid_pp(i)%mu_min = dflt_real fluid_pp(i)%mu_max = dflt_real fluid_pp(i)%mu_bulk = dflt_real + fluid_pp(i)%visc_model = visc_model_constant + fluid_pp(i)%visc_c = dflt_real + fluid_pp(i)%visc_d = dflt_real end do ! Subgrid bubble parameters diff --git a/src/simulation/m_arrhenius_viscosity.fpp b/src/simulation/m_arrhenius_viscosity.fpp new file mode 100644 index 0000000000..1c56b232ce --- /dev/null +++ b/src/simulation/m_arrhenius_viscosity.fpp @@ -0,0 +1,52 @@ +!> +!! @file m_arrhenius_viscosity.fpp +!! @brief Temperature-dependent (Arrhenius) shear viscosity mu = exp(C + D/T) and the mixture shear Reynolds number it gives. + +#:include 'macros.fpp' + +!> @brief Temperature-dependent (Arrhenius) shear viscosity mu = exp(C + D/T) and the mixture shear Reynolds number it gives. +module m_arrhenius_viscosity + + use m_derived_types !< Definitions of the derived types + use m_global_parameters !< sgm_eps, visc_models, visc_cs, visc_ds, fluid_inv_re + use m_constants, only: visc_exp_arg_cap + use m_eos, only: f_mixture_temperature + + implicit none + + private; public :: f_compute_arrhenius_shear_re + +contains + + !> Mixture shear Reynolds number 1/sum(alpha_i*mu_i) of one state at its thermal-equilibrium mixture temperature, the same + !! closure heat conduction uses. Arrhenius fluids take mu_i = exp(C_i + D_i/T); the others keep their constant 1/Re_i(1). + function f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) result(Re_shear) + + $:GPU_ROUTINE(function_name='f_compute_arrhenius_shear_re', parallelism='[seq]', cray_inline=True) + + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: alpha, alpha_rho + #:else + real(wp), dimension(num_fluids), intent(in) :: alpha, alpha_rho + #:endif + real(wp), intent(in) :: pres, gamma, pi_inf + real(wp) :: Re_shear + real(wp) :: T, inv_re + integer :: fl + + T = max(f_mixture_temperature(alpha_rho, pres, gamma, pi_inf), sgm_eps) + + inv_re = 0._wp + do fl = 1, num_fluids + if (visc_models(fl) == visc_model_arrhenius) then + inv_re = inv_re + alpha(fl)*exp(min(visc_cs(fl) + visc_ds(fl)/T, visc_exp_arg_cap)) + else + inv_re = inv_re + alpha(fl)*fluid_inv_re(fl) + end if + end do + + Re_shear = 1._wp/max(inv_re, sgm_eps) + + end function f_compute_arrhenius_shear_re + +end module m_arrhenius_viscosity diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 112ea52a23..01f5ae1f0d 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -16,6 +16,7 @@ module m_data_output use m_helper use m_helper_basic use m_sim_helpers + use m_arrhenius_viscosity use m_delay_file_access use m_ibm use m_boundary_common @@ -237,6 +238,8 @@ contains Re(1) = 1._wp/max(Re(1), sgm_eps) end if + if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) + call s_compute_stability_from_dt(vel, c, rho, Re, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) icfl_max_loc = max(icfl_max_loc, icfl) @@ -398,6 +401,8 @@ contains Re(1) = 1._wp/max(Re(1), sgm_eps) end if + if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) + call s_compute_stability_from_dt(vel, c, rho, Re, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) if (.not. f_approx_equal(icfl, icfl)) then diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 6157db87ef..a230460d6a 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -174,6 +174,15 @@ module m_global_parameters $:GPU_DECLARE(create='[any_non_newtonian, is_non_newtonian, hb_tau0, hb_K, hb_nn, hb_m_arr, hb_mu_min, hb_mu_max, fluid_inv_re]') + !> @name Temperature-dependent (Arrhenius) shear viscosity: per-fluid model ids and coefficients. + !> @{ + logical :: any_arrhenius_visc !< .true. if any fluid uses the Arrhenius mu(T) model + integer, allocatable, dimension(:) :: visc_models !< per-fluid shear viscosity model (visc_model_*) + real(wp), allocatable, dimension(:) :: visc_cs, visc_ds !< per-fluid Arrhenius coefficients C, D + !> @} + + $:GPU_DECLARE(create='[any_arrhenius_visc, visc_models, visc_cs, visc_ds]') + ! WENO averaging flag: use arithmetic mean or unaltered WENO-reconstructed cell-boundary values !> @{ real(wp) :: wa_flg @@ -401,6 +410,7 @@ contains shear_stress = .false. bulk_stress = .false. any_non_newtonian = .false. + any_arrhenius_visc = .false. num_igr_iters = dflt_num_igr_iters num_igr_warm_start_iters = dflt_num_igr_warm_start_iters alf_factor = dflt_alf_factor @@ -484,6 +494,9 @@ contains fluid_pp(i)%mu_min = dflt_real fluid_pp(i)%mu_max = dflt_real fluid_pp(i)%mu_bulk = dflt_real + fluid_pp(i)%visc_model = visc_model_constant + fluid_pp(i)%visc_c = dflt_real + fluid_pp(i)%visc_d = dflt_real end do ! Subgrid bubble parameters (bub_pp struct + scalar companions; scalar companions are @@ -889,6 +902,16 @@ contains end do $:GPU_UPDATE(device='[any_non_newtonian, is_non_newtonian, hb_tau0, hb_K, hb_nn, hb_m_arr, hb_mu_min, hb_mu_max, fluid_inv_re]') + ! Arrhenius mu(T) shear viscosity: gather per-fluid parameters into device arrays + @:ALLOCATE(visc_models(1:num_fluids), visc_cs(1:num_fluids), visc_ds(1:num_fluids)) + do i = 1, num_fluids + visc_models(i) = fluid_pp(i)%visc_model + visc_cs(i) = fluid_pp(i)%visc_c + visc_ds(i) = fluid_pp(i)%visc_d + end do + any_arrhenius_visc = any(visc_models == visc_model_arrhenius) + $:GPU_UPDATE(device='[any_arrhenius_visc, visc_models, visc_cs, visc_ds]') + if (bubbles_euler .and. qbmm .and. .not. polytropic) then allocate (MPI_IO_DATA%view(1:sys_size + 2*nb*nnode)) allocate (MPI_IO_DATA%var(1:sys_size + 2*nb*nnode)) @@ -1083,6 +1106,7 @@ contains @:DEALLOCATE(hb_tau0, hb_K, hb_nn, hb_m_arr) @:DEALLOCATE(hb_mu_min, hb_mu_max) @:DEALLOCATE(fluid_inv_re) + @:DEALLOCATE(visc_models, visc_cs, visc_ds) if (bubbles_euler) then @:DEALLOCATE(ptil) diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 5cddef70c3..0fe0e675df 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -18,6 +18,7 @@ module m_riemann_solver_hll use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights use m_thermochem_state, only: get_mixture_caloric_state use m_riemann_state + use m_arrhenius_viscosity implicit none @@ -200,6 +201,10 @@ contains if (viscous) then call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) call s_compute_interface_reynolds(alpha_R, Re_R, Re_size_loc1, Re_size_loc2) + if (any_arrhenius_visc) then + Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, alpha_rho_L, pres_L, gamma_L, pi_inf_L) + Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, alpha_rho_R, pres_R, gamma_R, pi_inf_R) + end if end if if (chemistry) then diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index 7bb7eec1a3..f6a3d814b0 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -21,6 +21,7 @@ module m_riemann_solver_hllc use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights use m_thermochem_state, only: get_mixture_caloric_state use m_riemann_state + use m_arrhenius_viscosity implicit none @@ -256,6 +257,10 @@ contains if (viscous) then call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) call s_compute_interface_reynolds(alpha_R, Re_R, Re_size_loc1, Re_size_loc2) + if (any_arrhenius_visc) then + Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, alpha_rho_L, pres_L, gamma_L, pi_inf_L) + Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, alpha_rho_R, pres_R, gamma_R, pi_inf_R) + end if end if call s_compute_energy(pres_L, alpha_rho_L, alpha_L, vel_L_rms, E_L) @@ -992,6 +997,10 @@ contains if (viscous) then call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) call s_compute_interface_reynolds(alpha_R, Re_R, Re_size_loc1, Re_size_loc2) + if (any_arrhenius_visc) then + Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, alpha_rho_L, pres_L, gamma_L, pi_inf_L) + Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, alpha_rho_R, pres_R, gamma_R, pi_inf_R) + end if end if if (chemistry) then diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index a918fec7cd..c2b4983242 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -25,6 +25,7 @@ module m_time_steppers use m_boundary_primitives, only: f_vel_ramp, bc_vel_ramp use m_helper use m_sim_helpers + use m_arrhenius_viscosity use m_fftw use m_nvtx use m_thermochem, only: num_species @@ -733,6 +734,8 @@ contains Re(1) = 1._wp/max(Re(1), sgm_eps) end if + if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) + call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, alpha, alpha_rho, j, k, l) ! min() drops a NaN candidate, so a cell with a NaN or non-positive dt would otherwise vanish from diff --git a/tests/6D29F444/golden-metadata.txt b/tests/6D29F444/golden-metadata.txt new file mode 100644 index 0000000000..89296923a5 --- /dev/null +++ b/tests/6D29F444/golden-metadata.txt @@ -0,0 +1,196 @@ +This file was created on 2026-10-23 13:24:15.130148. + +mfc.sh: + + Invocation: test --generate --only 6D29F444 -j 8 --no-debug -- -b mpirun + Lock: mpi=Yes & gpu=No & debug=No & reldebug=No & gcov=No & unified=No & single=No & mixed=No & fastmath=No + Git: 3ef5e1ee9bbecbfcdaef18cfb151fb304e438290 on pr/variable-viscosity-v2 (dirty) + +simulation: + + CMake Configuration: + + CMake v3.31.6 on ubuntu + + C : GNU v15.2.0 (/usr/bin/gcc) + Fortran : GNU v15.2.0 (/usr/bin/gfortran) + + PRE_PROCESS : OFF + SIMULATION : ON + POST_PROCESS : OFF + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /home/daveygravy/repos/MFC/.claude/worktrees/up-visc/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /usr/bin/gcc + CXX : /usr/bin/g++ + FC : /usr/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +post_process: + + CMake Configuration: + + CMake v3.31.6 on ubuntu + + C : GNU v15.2.0 (/usr/bin/gcc) + Fortran : GNU v15.2.0 (/usr/bin/gfortran) + + PRE_PROCESS : OFF + SIMULATION : OFF + POST_PROCESS : ON + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /home/daveygravy/repos/MFC/.claude/worktrees/up-visc/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /usr/bin/gcc + CXX : /usr/bin/g++ + FC : /usr/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +syscheck: + + CMake Configuration: + + CMake v3.31.6 on ubuntu + + C : GNU v15.2.0 (/usr/bin/gcc) + Fortran : GNU v15.2.0 (/usr/bin/gfortran) + + PRE_PROCESS : OFF + SIMULATION : OFF + POST_PROCESS : OFF + SYSCHECK : ON + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /home/daveygravy/repos/MFC/.claude/worktrees/up-visc/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /usr/bin/gcc + CXX : /usr/bin/g++ + FC : /usr/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +pre_process: + + CMake Configuration: + + CMake v3.31.6 on ubuntu + + C : GNU v15.2.0 (/usr/bin/gcc) + Fortran : GNU v15.2.0 (/usr/bin/gfortran) + + PRE_PROCESS : ON + SIMULATION : OFF + POST_PROCESS : OFF + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : OFF + + Fypp : /home/daveygravy/repos/MFC/.claude/worktrees/up-visc/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /usr/bin/gcc + CXX : /usr/bin/g++ + FC : /usr/bin/gfortran + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +CPU: + + CPU Info: + From lscpu + Architecture: x86_64 + CPU op-mode(s): 32-bit, 64-bit + Address sizes: 48 bits physical, 48 bits virtual + Byte Order: Little Endian + CPU(s): 256 + On-line CPU(s) list: 0-255 + Vendor ID: AuthenticAMD + Model name: AMD EPYC 7773X 64-Core Processor + CPU family: 25 + Model: 1 + Thread(s) per core: 2 + Core(s) per socket: 64 + Socket(s): 2 + Stepping: 2 + Frequency boost: enabled + CPU(s) scaling MHz: 83% + CPU max MHz: 3529.9060 + CPU min MHz: 401.4410 + BogoMIPS: 4400.14 + Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 invpcid cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin brs arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold v_vmsave_vmload vgif v_spec_ctrl umip pku ospke vaes vpclmulqdq rdpid overflow_recov succor smca debug_swap + Virtualization: AMD-V + L1d cache: 4 MiB (128 instances) + L1i cache: 4 MiB (128 instances) + L2 cache: 64 MiB (128 instances) + L3 cache: 1.5 GiB (16 instances) + NUMA node(s): 2 + NUMA node0 CPU(s): 0-63,128-191 + NUMA node1 CPU(s): 64-127,192-255 + Vulnerability Gather data sampling: Not affected + Vulnerability Ghostwrite: Not affected + Vulnerability Indirect target selection: Not affected + Vulnerability Itlb multihit: Not affected + Vulnerability L1tf: Not affected + Vulnerability Mds: Not affected + Vulnerability Meltdown: Not affected + Vulnerability Mmio stale data: Not affected + Vulnerability Old microcode: Not affected + Vulnerability Reg file data sampling: Not affected + Vulnerability Retbleed: Not affected + Vulnerability Spec rstack overflow: Vulnerable + Vulnerability Spec store bypass: Vulnerable + Vulnerability Spectre v1: Vulnerable: __user pointer sanitization and usercopy barriers only; no swapgs barriers + Vulnerability Spectre v2: Vulnerable; IBPB: disabled; STIBP: disabled; PBRSB-eIBRS: Not affected; BHI: Not affected + Vulnerability Srbds: Not affected + Vulnerability Tsa: Vulnerable + Vulnerability Tsx async abort: Not affected + Vulnerability Vmscape: Vulnerable + diff --git a/tests/6D29F444/golden.txt b/tests/6D29F444/golden.txt new file mode 100644 index 0000000000..8c511ea165 --- /dev/null +++ b/tests/6D29F444/golden.txt @@ -0,0 +1,16 @@ +D/cons.1.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 +D/cons.1.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000002 1.00000000000031 1.0000000000044 1.00000000005348 1.00000000054761 1.00000000462965 1.00000002487515 0.50000011143128 0.50000000676529 0.50000000138136 0.50000000026174 0.50000000004278 0.50000000000608 0.50000000000076 0.50000000000008 0.50000000000001 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.49999997215963 0.12500002784037 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 +D/cons.2.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +D/cons.2.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999987 0.99999999999675 0.99999999992663 0.999999998492 0.99999997210299 0.99999954110118 0.99999338962525 0.99991825392086 0.99915497202562 0.99296638971589 0.95541178306847 0.80206622717864 0.17094858929074 0.05878528809715 0.01634729542607 0.00362608218487 0.00066367310572 0.00010305836376 1.387243299e-05 1.64612646e-06 1.7452035e-07 1.67127e-08 1.45878e-09 1.1694e-10 8.67e-12 6e-13 4e-14 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2e-14 1.7e-13 1.45e-12 1.115e-11 7.511e-11 4.3423e-10 2.10735e-09 8.35264e-09 2.619791e-08 6.97319e-09 5.12267e-09 3.65936e-09 2.53878e-09 1.70876e-09 1.11472e-09 7.0427e-10 4.3065e-10 2.5475e-10 1.4574e-10 8.062e-11 4.311e-11 2.229e-11 1.115e-11 5.39e-12 2.52e-12 1.14e-12 5e-13 2.1e-13 9e-14 3e-14 1e-14 1e-14 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 +D/cons.3.00.000000.dat 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 +D/cons.3.00.000050.dat 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 2.99999999999987 2.99999999999676 2.99999999992663 2.999999998492 2.999999972103 2.99999954110147 2.9999933896742 2.99991826115835 2.99915573596077 2.99301949016276 2.9576393338427 2.84619608664533 1.44371345204318 1.2596590195755 1.25067154824429 1.2500335684872 1.2500011723089 1.25000002969998 1.25000000056833 1.25000000000867 1.25000000000013 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.24999992195002 0.25000007804998 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 +D/cons.4.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 +D/cons.4.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 +D/prim.1.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 +D/prim.1.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000002 1.00000000000031 1.0000000000044 1.00000000005348 1.00000000054761 1.00000000462965 1.00000002487515 0.50000011143128 0.50000000676529 0.50000000138136 0.50000000026174 0.50000000004278 0.50000000000608 0.50000000000076 0.50000000000008 0.50000000000001 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.49999997215963 0.12500002784037 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 +D/prim.2.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +D/prim.2.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.99999999999987 0.99999999999675 0.99999999992663 0.999999998492 0.99999997210299 0.99999954110116 0.99999338962494 0.99991825391646 0.99915497197218 0.99296638917213 0.95541177864525 0.80206620722713 0.34189710238542 0.11757057460351 0.03269459076182 0.00725216436594 0.00132734621132 0.00020611672751 2.774486597e-05 3.29225291e-06 3.490407e-07 3.34254e-08 2.91756e-09 2.3388e-10 1.733e-11 1.19e-12 8e-14 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3e-14 3.3e-13 2.89e-12 2.231e-11 1.5022e-10 8.6846e-10 4.2147e-09 1.670528e-08 5.239582e-08 5.578554e-08 4.098135e-08 2.927487e-08 2.03102e-08 1.367007e-08 8.91776e-09 5.63416e-09 3.44522e-09 2.03804e-09 1.16593e-09 6.4493e-10 3.4491e-10 1.7834e-10 8.917e-11 4.313e-11 2.018e-11 9.14e-12 4.01e-12 1.7e-12 7e-13 2.8e-13 1.1e-13 4e-14 1e-14 1e-14 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 +D/prim.3.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 +D/prim.3.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00000000000008 1.0000000000109 1.00000000155939 1.00000016277028 1.00001134595201 1.00049339933707 1.0098163913025 0.56579601535019 0.50248132381024 0.50016172567091 0.50000816800608 0.50000029273876 0.50000000763158 0.50000000015035 0.50000000000239 0.50000000000004 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.49999996878001 0.10000003121999 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 +D/prim.4.00.000000.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 +D/prim.4.00.000050.dat 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 \ No newline at end of file diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 476009fbd0..0ae665fc46 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -222,6 +222,17 @@ "does not have. Not supported with igr or chemistry (which carries its own mixture-averaged conduction)." ), }, + "check_visc_model": { + "title": "Temperature-Dependent Viscosity", + "category": "Numerical Schemes", + "math": r"\mu_k(T) = \exp\!\left(C_{\mu,k} + D_{\mu,k}/T\right)", + "explanation": ( + "fluid_pp(i)%visc_model = 1 (Arrhenius) requires fluid_pp(i)%visc_c and a nonzero fluid_pp(i)%visc_d; " + "visc_c and visc_d are rejected on other fluids. Requires viscous = T, model_eqns = 2 or 3, " + "riemann_solver 1 (HLL) or 2 (HLLC), fluid_pp(i)%cv > 0 on every fluid, and the stiffened-gas or ideal-gas " + "equation of state. Not supported with non_newtonian, cyl_coord, bubbles_euler, ib, igr, or chemistry." + ), + }, # Feature Compatibility "check_mhd": { "title": "Magnetohydrodynamics (MHD)", @@ -1560,6 +1571,62 @@ def check_heat_conduction(self): "heat conduction is not supported with chemistry: the reacting path already carries mixture-averaged conduction through chem_params%diffusion", ) + def check_visc_model(self): + """Checks constraints on the temperature-dependent (Arrhenius) shear viscosity model (fluid_pp(i)%visc_model)""" + num_fluids = self.get("num_fluids") + # If num_fluids is not set, check at least fluid 1 (for model_eqns=1) + if num_fluids is None: + num_fluids = 1 + arrhenius = CONSTRAINTS["fluid_pp(1)%visc_model"]["names"]["arrhenius"] + eos_names = CONSTRAINTS["fluid_pp(1)%eos"]["names"] + supported_eos = {eos_names["stiffened_gas"], eos_names["ideal_gas"]} + + any_arrhenius = False + for i in range(1, num_fluids + 1): + if self.get(f"fluid_pp({i})%visc_model") != arrhenius: + for coeff in ["visc_c", "visc_d"]: + self.prohibit( + self.get(f"fluid_pp({i})%{coeff}") is not None, + f"fluid_pp({i})%{coeff} is set, but fluid_pp({i})%visc_model is not the Arrhenius model", + ) + continue + + any_arrhenius = True + visc_d = self.get(f"fluid_pp({i})%visc_d") + self.prohibit(self.get(f"fluid_pp({i})%visc_c") is None, f"fluid_pp({i})%visc_model = Arrhenius requires fluid_pp({i})%visc_c to be set") + self.prohibit(visc_d is None or visc_d == 0, f"fluid_pp({i})%visc_model = Arrhenius requires a nonzero fluid_pp({i})%visc_d (D = 0 is a constant viscosity)") + + if not any_arrhenius: + return + + model_eqns = self.get("model_eqns") + riemann_solver = self.get("riemann_solver") + self.prohibit(self.get("viscous", "F") != "T", "the Arrhenius viscosity model requires viscous = T") + self.prohibit(model_eqns not in (2, 3), "the Arrhenius viscosity model requires model_eqns = 2 (5-equation) or 3 (6-equation)") + self.prohibit( + riemann_solver is not None and riemann_solver not in (1, 2), + "the Arrhenius viscosity model requires riemann_solver 1 (HLL) or 2 (HLLC)", + ) + for i in range(1, num_fluids + 1): + cv = self.get(f"fluid_pp({i})%cv") + # The mixture temperature divides by sum(alpha_rho_i*cv_i*n_i); a fluid with cv = 0 would drop out of it + self.prohibit(cv is None or cv <= 0, f"fluid_pp({i})%cv must be positive on every fluid when the Arrhenius viscosity model is used") + eos = self.get(f"fluid_pp({i})%eos") + self.prohibit( + (eos if eos is not None else eos_names["stiffened_gas"]) not in supported_eos, + "the Arrhenius viscosity model supports only the stiffened-gas and ideal-gas equations of state", + ) + self.prohibit( + self.get(f"fluid_pp({i})%non_newtonian", "F") == "T", + "the Arrhenius viscosity model cannot be combined with non_newtonian: the Herschel-Bulkley stress does not see mu(T)", + ) + # These paths evaluate the shear stress from the constant fluid_pp(i)%Re(1) rather than the face mu(T) + self.prohibit(self.get("cyl_coord", "F") == "T", "the Arrhenius viscosity model is not supported with cyl_coord") + self.prohibit(self.get("bubbles_euler", "F") == "T", "the Arrhenius viscosity model is not supported with bubbles_euler") + self.prohibit(self.get("ib", "F") == "T", "the Arrhenius viscosity model is not supported with ib") + self.prohibit(self.get("igr", "F") == "T", "the Arrhenius viscosity model is not supported with igr") + self.prohibit(self.get("chemistry", "F") == "T", "the Arrhenius viscosity model is not supported with chemistry") + def check_non_newtonian(self): """Checks constraints on non-Newtonian (Herschel-Bulkley) parameters (simulation)""" num_fluids = self.get("num_fluids") @@ -3059,6 +3126,7 @@ def validate_simulation(self): self.check_synthetic_turbulence() self.check_viscosity() self.check_heat_conduction() + self.check_visc_model() self.check_non_newtonian() self.check_mhd_simulation() self.check_igr_simulation() diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index ba5613db24..45d71d1600 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -955,6 +955,15 @@ def _load(): _r(f"{px}mu_min", REAL, {"viscosity"}, math=r"\mu_{\min,k}") _r(f"{px}mu_max", REAL, {"viscosity"}, math=r"\mu_{\max,k}") _r(f"{px}mu_bulk", REAL, {"viscosity"}, math=r"\mu_{\mathrm{bulk},k}") + CONSTRAINTS[f"fluid_pp({f})%visc_model"] = { + "choices": [0, 1], + "value_labels": {0: "constant", 1: "Arrhenius mu(T)"}, + "names": {"constant": 0, "arrhenius": 1}, + "fortran_prefix": "visc_model", + } + _r(f"{px}visc_model", INT, {"viscosity"}) + _r(f"{px}visc_c", REAL, {"viscosity"}, math=r"C_{\mu,k}") + _r(f"{px}visc_d", REAL, {"viscosity"}, math=r"D_{\mu,k}") # bub_pp (bubble properties) for a, sym in [ diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 7dfb46f3ea..588e036aff 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -363,6 +363,9 @@ (r"fluid_pp\((\d+)\)%mu_min", "Lower viscosity clamp for fluid {0}"), (r"fluid_pp\((\d+)\)%mu_max", "Upper viscosity clamp for fluid {0}"), (r"fluid_pp\((\d+)\)%mu_bulk", "Bulk viscosity (non-Newtonian) for fluid {0}"), + (r"fluid_pp\((\d+)\)%visc_model", "Shear viscosity model for fluid {0} (0 constant, 1 Arrhenius mu(T))"), + (r"fluid_pp\((\d+)\)%visc_c", "Arrhenius viscosity coefficient C for fluid {0}"), + (r"fluid_pp\((\d+)\)%visc_d", "Arrhenius viscosity coefficient D for fluid {0}"), # patch_ib patterns (r"patch_ib\((\d+)\)%geometry", "Geometry type for immersed boundary {0}"), (r"patch_ib\((\d+)\)%x_centroid", "X-coordinate of centroid for IB patch {0}"), diff --git a/toolchain/mfc/params_tests/test_eos_families.py b/toolchain/mfc/params_tests/test_eos_families.py index d6d3514b75..95e370208b 100644 --- a/toolchain/mfc/params_tests/test_eos_families.py +++ b/toolchain/mfc/params_tests/test_eos_families.py @@ -65,6 +65,9 @@ "fluid_pp(1)%vinet_k0p", "fluid_pp(1)%vinet_rho0", "fluid_pp(1)%vinet_t0", + "fluid_pp(1)%visc_c", + "fluid_pp(1)%visc_d", + "fluid_pp(1)%visc_model", ] diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index d169be52bc..4ceb53b887 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -883,6 +883,10 @@ def alter_num_fluids(dimInfo): cases.append(define_case_d(stack, "weno_avg", {"weno_avg": "T"})) stack.pop() + if len(dimInfo[0]) == 1: + # Arrhenius mu(T): the patches start at T = 2.5, 2.5, 2.0, so mu spans e^9 to e^10 around the base 1/Re(1) = 1e4 + cases.append(define_case_d(stack, "visc_model=1", {"fluid_pp(1)%visc_model": 1, "fluid_pp(1)%visc_c": 5.0, "fluid_pp(1)%visc_d": 10.0, "fluid_pp(1)%cv": 1.0})) + stack.pop() if len(dimInfo[0]) <= 2: diff --git a/toolchain/mfc/test_case_validator.py b/toolchain/mfc/test_case_validator.py index c0170452d9..2a60167b98 100644 --- a/toolchain/mfc/test_case_validator.py +++ b/toolchain/mfc/test_case_validator.py @@ -704,3 +704,40 @@ def test_rejects_isothermal_wall_without_any_heat_path(self): {**BASE, "bc_x%beg": -16, "bc_x%end": -16, "bc_x%isothermal_in": "T", "bc_x%Twall_in": 300.0}, "requires a heat-conduction path", ) + + +class TestArrheniusViscosity(ConstraintTestCase): + """fluid_pp(i)%visc_model = 1 (mu(T) = exp(visc_c + visc_d/T)) input rules. The model needs both + coefficients, a nonzero D, and the mixture temperature (cv > 0, a supported EOS); it is only wired + into the HLL/HLLC Riemann states and the viscous dt, so the paths that bypass them are refused.""" + + GOOD = {**BASE, "viscous": "T", "fluid_pp(1)%Re(1)": 100.0, "fluid_pp(1)%cv": 1.0, "fluid_pp(1)%visc_model": 1, "fluid_pp(1)%visc_c": 0.0, "fluid_pp(1)%visc_d": 1.0} + + def test_accepts_valid_configuration(self): + self.assertAccepts(self.GOOD) + + def test_rejects_missing_coefficients(self): + self.assertRejects({**self.GOOD, "fluid_pp(1)%visc_c": None}, "requires fluid_pp(1)%visc_c to be set") + self.assertRejects({**self.GOOD, "fluid_pp(1)%visc_d": None}, "requires a nonzero fluid_pp(1)%visc_d") + + def test_rejects_zero_d(self): + """D = 0 is a constant viscosity exp(C); the constant model already says that.""" + self.assertRejects({**self.GOOD, "fluid_pp(1)%visc_d": 0.0}, "requires a nonzero fluid_pp(1)%visc_d") + + def test_rejects_coefficients_without_the_model(self): + self.assertRejects({**self.GOOD, "fluid_pp(1)%visc_model": 0}, "fluid_pp(1)%visc_c is set, but fluid_pp(1)%visc_model is not the Arrhenius model") + + def test_rejects_missing_cv(self): + self.assertRejects({**self.GOOD, "fluid_pp(1)%cv": None}, "fluid_pp(1)%cv must be positive on every fluid") + + def test_rejects_without_viscous(self): + self.assertRejects({**self.GOOD, "viscous": "F"}, "requires viscous = T") + + def test_rejects_other_riemann_solvers(self): + self.assertRejects({**self.GOOD, "riemann_solver": 5}, "requires riemann_solver 1 (HLL) or 2 (HLLC)") + + def test_rejects_non_newtonian(self): + self.assertRejects({**self.GOOD, "fluid_pp(1)%non_newtonian": "T"}, "cannot be combined with non_newtonian") + + def test_rejects_cyl_coord(self): + self.assertRejects({**self.GOOD, "cyl_coord": "T"}, "not supported with cyl_coord") From 0191d971476b403420cb8851a3b19998b0ce8f79 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Fri, 23 Oct 2026 14:43:30 -0400 Subject: [PATCH 2/4] Add 2D Couette variable-viscosity validation example Plane Couette flow between a cold fixed wall and a hot sliding wall with Arrhenius mu(T), validated against the exact steady solution of the coupled momentum/energy BVP. Uses upstream k_therm conduction with isothermal no-slip walls. Observed L2 velocity order 2.00 over Ny = 33, 65, 97. Skipped in the example test sweep: its analytic density IC needs a case-specific pre_process build; mu(T) is covered by golden 6D29F444. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VtsGVvCfsDJz5dmtSk9dKt --- .../2D_Couette_Variable_Viscosity/.gitignore | 6 + .../2D_Couette_Variable_Viscosity/README.md | 155 +++++++++++++ .../2D_Couette_Variable_Viscosity/case.py | 145 ++++++++++++ .../couette_config.py | 101 +++++++++ .../reference.py | 97 ++++++++ .../run_suite.sh | 56 +++++ .../2D_Couette_Variable_Viscosity/validate.py | 208 ++++++++++++++++++ toolchain/mfc/test/cases.py | 1 + 8 files changed, 769 insertions(+) create mode 100644 examples/2D_Couette_Variable_Viscosity/.gitignore create mode 100644 examples/2D_Couette_Variable_Viscosity/README.md create mode 100644 examples/2D_Couette_Variable_Viscosity/case.py create mode 100644 examples/2D_Couette_Variable_Viscosity/couette_config.py create mode 100644 examples/2D_Couette_Variable_Viscosity/reference.py create mode 100755 examples/2D_Couette_Variable_Viscosity/run_suite.sh create mode 100644 examples/2D_Couette_Variable_Viscosity/validate.py diff --git a/examples/2D_Couette_Variable_Viscosity/.gitignore b/examples/2D_Couette_Variable_Viscosity/.gitignore new file mode 100644 index 0000000000..ac6874d94d --- /dev/null +++ b/examples/2D_Couette_Variable_Viscosity/.gitignore @@ -0,0 +1,6 @@ +# MFC run outputs and the per-grid run copies live here; regenerate with run_suite.sh. +runs/ +*.log +# validate.py output +figures/ +summary.json diff --git a/examples/2D_Couette_Variable_Viscosity/README.md b/examples/2D_Couette_Variable_Viscosity/README.md new file mode 100644 index 0000000000..a5e91ebaec --- /dev/null +++ b/examples/2D_Couette_Variable_Viscosity/README.md @@ -0,0 +1,155 @@ +# 2D Couette flow with temperature-dependent viscosity + +A unit-level validation of MFC's Arrhenius viscosity model +`fluid_pp%visc_model = 1`, `mu(T) = exp(C + D/T)`, against an **exact** steady +solution. This is the analytical gold standard for the variable-viscosity +feature: it pins down the spatial order of accuracy on a problem with a known +closed-form answer, so a sign error or discretization bug in the variable-`mu` +viscous flux shows up directly. + +## The physics, and why it tests mu(T) + +A single fluid fills the gap `0 <= y <= H` between two no-slip walls: + +- bottom wall (`y = 0`): fixed, held at temperature `T0 = 300 K` (cold), +- top wall (`y = H`): sliding at speed `U`, held at temperature `T1 = 400 K` (hot), +- periodic in the streamwise (`x`) direction. + +Conduction sustains a temperature gradient across the gap, and because +`mu(T) = exp(C + D/T)` with `D > 0`, the viscosity **falls as temperature rises** +(the liquid/Arrhenius law). The hot fluid near the moving wall is therefore +*thinner* than the cold fluid near the fixed wall. + +In steady state the shear stress `mu(T) du/dy` is uniform across the gap (no +pressure gradient, no body force). Uniform stress with a non-uniform `mu` forces +`du/dy` to vary inversely with `mu`: the velocity climbs slowly through the cold, +viscous lower half and steeply through the hot, thin upper half. **The velocity +profile is curved.** A constant-viscosity solver gives a straight line +`u = U y/H` with `u(H/2) = U/2`; here the exact mid-gap velocity is +`u(H/2) ≈ 0.352 U`. That ~30% deficit *is* the `mu(T)` signal, and it must match +the exact profile to within the discretization error. + +## The exact solution (`reference.py`) + +The incompressible steady state solves the coupled two-point BVP + +``` +momentum: d/dy [ mu(T) du/dy ] = 0 +energy: d/dy [ k dT/dy ] + mu(T) (du/dy)^2 = 0 +``` + +with `u(0)=0, u(H)=U, T(0)=T0, T(H)=T1`. `reference.py` solves it with +`scipy.solve_bvp` to ~1e-10. The viscous-heating term `mu(T)(du/dy)^2` couples +the energy equation back to the momentum equation, so the reference also exercises +the energy side of `mu(T)`. At this case's Brinkman number (`Br ≈ 0.014`) the +heating bump is small (~0.17 K above the linear profile), and the closed-form +`Br -> 0` limit (linear `T`, `u = U * integral_0^y dy'/mu / integral_0^H dy'/mu`) +agrees with the coupled solver to ~0.02 % of `U` -- an independent check on the +solver. + +## Compressible-solver setup + +MFC is a compressible code; temperature is not stored but follows from the +stiffened-gas EOS `T = (p + p_inf)/((gam-1) rho cv)`. Two consequences shape the +deck (`case.py`, constants in `couette_config.py`): + +- **Soft sound speed.** A deliberately small `cv` puts the sound speed near + `c ≈ 10` so the flow sits at `Ma ≈ 0.1`. Low Mach keeps the incompressible + reference faithful and the acoustic time step affordable. +- **Temperature via a density profile.** The target linear `T(y)` is imposed by + initializing a `y`-varying density `rho(y) = (p0+p_inf)/((gam-1) cv T(y))` at + uniform pressure -- a valid mechanical equilibrium (no gravity). +- **Modes.** The deck uses `model_eqns = 3` (6-equation) and `riemann_solver = 2` + (HLLC); `visc_model = 1` also runs with `model_eqns = 2` and HLL. MFC evaluates + `mu(T)` at the mixture temperature of the reconstructed left and right states at + each cell face, so the face viscosity follows the reconstructed temperature. + Fourier conduction (`fluid_pp(1)%k_therm > 0`) with isothermal no-slip walls + (`bc_y = -16`, `isothermal_in/out`) holds the walls at `T0` and `T1`. +- **Compact viscous gradient (`weno_Re_flux = F`, `weno_avg = F`).** `weno_Re_flux` + is designed for *discontinuous* viscosity at material interfaces. For a single + fluid with smooth `mu(T)` it leaves the `2*dy` (odd-even) mode of the + wall-normal velocity undamped -- and because this flow has no advection + (`v ~ 0`), the HLLC upwinding has nothing to damp either -- so `u(y)` grows a + grid-scale checkerboard (~5 % of `U` at `Ny = 33`). The compact central + gradient damps it; the resolved profile then matches the exact solution to the + discretization error. (The two flags must be set together: `weno_avg = F` with + `weno_Re_flux = T` disables the viscous coupling entirely.) + +The flow starts from rest (`u = 0`) and develops toward the curved steady profile, +so the match is *not* baked into the initial condition. + +Dimensionless parameters: `Re ≈ 20`, `Ma ≈ 0.1`, `Pr = 1`, `Br ≈ 0.014`, +viscosity contrast `mu(T0)/mu(T1) ≈ 3.5`. + +## Running it + +```bash +# 1. Run the grid-refinement sweep. The first grid is a normal building run +# (the analytic density IC compiles into pre_process, so the build is +# case-specific); the remaining grids then reuse that binary (--no-build) +# and run concurrently, each in its own runs/n/ on one rank. +# Default grids: Ny = 33, 65, 97. +cd examples/2D_Couette_Variable_Viscosity +./run_suite.sh # or: ./run_suite.sh 32 64 96 128 + +# 2. Compare to the exact solution; writes summary.json + figures/ +python3 validate.py +``` + +A single grid can be run directly, choosing the wall-normal resolution with an +environment variable: + +```bash +COUETTE_N=64 ./mfc.sh run examples/2D_Couette_Variable_Viscosity/case.py -n 1 +``` + +## Files + +| file | role | +|------|------| +| `couette_config.py` | single source of truth for all physical constants | +| `reference.py` | exact coupled-BVP solution (`scipy.solve_bvp`) | +| `case.py` | MFC input deck, parametrized by `COUETTE_N` | +| `run_suite.sh` | runs the grid sweep into `runs/n/` | +| `validate.py` | compares MFC to the exact solution; figures + `summary.json` | + +## Results + +The MFC steady velocity matches the exact `mu(T)` Couette profile and converges at +second order; temperature is exact to ~1e-6 and the runs are fully steady +(`unsteadiness ~ 1e-9`). The viscosity contrast bends `u(H/2)` to **0.352** versus +the constant-viscosity **0.500** -- the signal a constant-`mu` code cannot produce. + +| Ny | L2(u)/U | Linf(u)/U | L2(T)/dT | +|----|----------|-----------|----------| +| 33 | 1.33e-4 | 1.83e-4 | 1.58e-6 | +| 65 | 3.41e-5 | 4.85e-5 | 3.93e-7 | +| 97 | 1.53e-5 | 2.20e-5 | 1.75e-7 | + +Observed spatial order (L2 velocity) = **2.00**. `validate.py` writes to `figures/`: profile +overlay (`couette_profiles.png`), L2 convergence (`couette_convergence.png`), and +the spatial relative error (`couette_error.png`) -- the latter shows the +pointwise `(u - u_exact)/U` shrinking with refinement and collapsing under an +`Ny^2` rescaling, i.e. the second-order accuracy holds everywhere in the gap, not +just in the L2 norm, and the full data to `summary.json`. + +## Scope and the companion benchmark + +This case validates the **exponential/Arrhenius (liquid) law** that MFC actually +implements -- the right gold standard for temperature-dependent droplet +viscosity. It is a unit test: `mu(T)` is isolated on an exact answer, free of the +EOS/buoyancy coupling that entangles system-level benchmarks. + +The natural *system-level* companion is the low-Mach differentially-heated cavity +(Le Quéré / Vierendeels), where boundary-layer asymmetry is driven by the +viscosity contrast and validated against tabulated wall Nusselt numbers. That +benchmark uses **Sutherland's law (a gas law)**, which MFC does not yet +implement, so it is a separate follow-up (add the Sutherland model, then the +cavity case). + +### Cost note + +The time step is explicit-stability limited: acoustic (`~dy`) on coarse grids and +**viscous (`~dy^2`)** once `Ny` exceeds ~60. The default sweep stays in the cheap +range; `Ny = 129` is viscous-CFL limited (~221k steps, hours on a single rank) and +is optional. diff --git a/examples/2D_Couette_Variable_Viscosity/case.py b/examples/2D_Couette_Variable_Viscosity/case.py new file mode 100644 index 0000000000..6878a978c2 --- /dev/null +++ b/examples/2D_Couette_Variable_Viscosity/case.py @@ -0,0 +1,145 @@ +#!/usr/bin/env python3 +"""Variable-viscosity planar Couette flow -- MFC input deck. + +Validates the temperature-dependent (Arrhenius) viscosity feature +fluid_pp%visc_model = 1, mu(T) = exp(C + D/T), against the exact steady solution +in reference.py. See README.md for the full description. + +Setup (all constants from couette_config.py): + - 2D channel, periodic in x, walls top and bottom in y. + - Bottom wall (y=0): no-slip, fixed, isothermal at T0 (cold -> viscous). + - Top wall (y=H): no-slip, sliding at U, isothermal at T1 (hot -> thin). + - Single fluid; model_eqns = 3 + riemann_solver = 2 (HLLC). + - viscous on, and Fourier conduction (fluid_pp(1)%k_therm > 0) pins the wall + temperatures and sustains the temperature gradient that drives mu(T). + - Initial condition: quiescent (u=0), uniform pressure, with a y-varying + density that encodes the linear conduction temperature profile via the EOS. + The flow then develops from rest toward the curved steady Couette profile. + +The wall-normal grid n is set by the COUETTE_N environment variable so the +validation harness can run a spatial-refinement sweep on a single deck. +""" + +import json +import os + +import couette_config as cfg + +# Wall-normal resolution (refined by the convergence sweep); streamwise grid is +# fixed and small because there is no x-gradient. +# WENO5 needs >= 25 cells per dimension, so the streamwise (x) grid is pinned at +# that minimum (there is no x-gradient) and the sweep keeps n >= 32. +n_cells = int(os.environ.get("COUETTE_N", "64")) +m_cells = int(os.environ.get("COUETTE_M", "24")) + +# Square cells: size the (periodic, gradient-free) x-extent so dx == dy on every +# grid. This keeps the acoustic CFL set by dy alone and the time step ~ 1/n. +dy = cfg.H / (n_cells + 1) +Lx = dy * (m_cells + 1) +dx = Lx / (m_cells + 1) + +# Time step: minimum of acoustic, viscous, and thermal explicit-stability limits. +c_max = (cfg.gam * (cfg.gam - 1.0) * cfg.cv * cfg.T1) ** 0.5 # sound speed, hot wall +dt_acoustic = 0.4 * min(dx, dy) / c_max +dt_viscous = 0.4 * dy**2 / (2.0 * cfg.nu_max) +alpha_max = cfg.k_therm / (cfg.rho_of_T(cfg.T0) * cfg.cp) # thermal diffusivity, cold wall +dt_thermal = 0.4 * dy**2 / (2.0 * alpha_max) +dt = min(dt_acoustic, dt_viscous, dt_thermal) + +t_step_stop = int(round(cfg.t_end / dt)) +if "COUETTE_NSTEPS" in os.environ: # short override for timing probes + t_step_stop = int(os.environ["COUETTE_NSTEPS"]) +t_step_save = max(1, t_step_stop // 40) # ~40 snapshots to confirm steadiness + +# Linear conduction temperature profile T(y) = T0 + (T1-T0)*y/H imposed through +# the density proxy rho(y) = (p0+p_inf)/((gam-1)*cv*T(y)). Written with plain +# decimals (no 'e' literals -- MFC would expand a bare 'e' to Euler's number). +gm1_cv = (cfg.gam - 1.0) * cfg.cv +slope = (cfg.T1 - cfg.T0) / cfg.H +rho_expr = f"({cfg.p0 + cfg.p_inf})/(({gm1_cv})*({cfg.T0} + ({slope})*y))" + +case = { + # -- Run info -- + "run_time_info": "F", + # -- Domain (periodic x, walls in y) -- + "x_domain%beg": 0.0, + "x_domain%end": Lx, + "y_domain%beg": 0.0, + "y_domain%end": cfg.H, + "m": m_cells, + "n": n_cells, + "p": 0, + "cyl_coord": "F", + # -- Time stepping -- + "dt": dt, + "t_step_start": 0, + "t_step_stop": t_step_stop, + "t_step_save": t_step_save, + # -- Model: 6-equation + HLLC (visc_model = 1 also runs with model_eqns = 2 or HLL) -- + "model_eqns": 3, + "num_fluids": 1, + "alt_soundspeed": "F", + "mpp_lim": "F", + "mixture_err": "T", + "time_stepper": 3, + # -- Spatial scheme -- + "weno_order": 5, + "weno_eps": 1e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + # Compact (non-WENO) viscous gradient. weno_Re_flux is meant for discontinuous + # viscosity at material interfaces; for a single fluid with smooth mu(T) it + # leaves the 2*dy (odd-even) mode of the wall-normal velocity undamped -- in + # this advection-free flow there is no convective upwinding to remove it -- so + # u(y) develops a grid-scale checkerboard. The compact central gradient damps + # it and the profile matches the exact solution. (Both must be set together: + # weno_avg=F with weno_Re_flux=T disables the viscous coupling entirely.) + "weno_Re_flux": os.environ.get("COUETTE_WENO_RE", "F"), + "weno_avg": os.environ.get("COUETTE_WENO_AVG", "F"), + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + # -- Boundary conditions: periodic x; no-slip isothermal walls in y -- + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -16, # no-slip wall, bottom (fixed) + "bc_y%end": -16, # no-slip wall, top (sliding) + "bc_y%ve1": cfg.U, # top-wall x-velocity (the Couette drive) + "bc_y%isothermal_in": "T", + "bc_y%Twall_in": cfg.T0, + "bc_y%isothermal_out": "T", + "bc_y%Twall_out": cfg.T1, + # -- Output -- + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "T", + # -- Physics -- + "viscous": "T", + "chemistry": "F", + # -- Fluid: ideal gas with Arrhenius mu(T) = exp(C + D/T) -- + "fluid_pp(1)%gamma": 1.0 / (cfg.gam - 1.0), + "fluid_pp(1)%pi_inf": cfg.gam * cfg.p_inf / (cfg.gam - 1.0), + "fluid_pp(1)%cv": cfg.cv, + "fluid_pp(1)%Re(1)": 1.0 / cfg.mu_ref, # nominal Re; visc_model=1 overrides with mu(T) + "fluid_pp(1)%k_therm": cfg.k_therm, + "fluid_pp(1)%visc_model": 1, + "fluid_pp(1)%visc_c": cfg.C, + "fluid_pp(1)%visc_d": cfg.D, + # -- Initial condition: quiescent, uniform p, linear-T density profile -- + "num_patches": 1, + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": Lx / 2, + "patch_icpp(1)%y_centroid": cfg.H / 2, + "patch_icpp(1)%length_x": Lx, + "patch_icpp(1)%length_y": cfg.H, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": cfg.p0, + "patch_icpp(1)%alpha_rho(1)": rho_expr, + "patch_icpp(1)%alpha(1)": 1.0, +} + +if __name__ == "__main__": + print(json.dumps(case)) diff --git a/examples/2D_Couette_Variable_Viscosity/couette_config.py b/examples/2D_Couette_Variable_Viscosity/couette_config.py new file mode 100644 index 0000000000..6ff37af8f9 --- /dev/null +++ b/examples/2D_Couette_Variable_Viscosity/couette_config.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python3 +"""Single source of truth for the variable-viscosity Couette validation. + +case.py (the MFC input), reference.py (the exact solution), and validate.py +(the comparison) all import these constants so they can never drift apart. + +Physics: planar Couette flow of a single fluid between two no-slip walls a +distance H apart. The bottom wall (y = 0) is held fixed at temperature T0; the +top wall (y = H) slides at speed U and is held at temperature T1. The fluid has +a temperature-dependent (Arrhenius) dynamic viscosity + + mu(T) = exp(C + D / T) with D > 0 -> mu falls as T rises + +so the hot side near the top wall is thinner than the cold side near the +bottom. In steady state the shear stress mu(T) du/dy is uniform across the gap, +which forces du/dy to vary inversely with mu -- the velocity profile is *curved* +rather than the straight line a constant-viscosity code would give. That +curvature is the signature of mu(T), and it is what this case validates against +an exact solution. + +MFC is a compressible solver, so we deliberately soften the sound speed (small +cv) to keep the flow at low Mach number; that makes the incompressible Couette +solution in reference.py a faithful reference and keeps the acoustic time step +affordable. Temperature is not a stored field: it follows from the stiffened-gas +EOS T = (p + p_inf)/((gam - 1) rho cv), so a target temperature profile is +imposed through a density profile rho(T) at uniform pressure. +""" + +import numpy as np + +# Gas / EOS +gam = 1.4 # ratio of specific heats +p_inf = 0.0 # ideal gas (no stiffening needed: the soft cv already softens c) + +# Geometry and drive (the streamwise x-extent is set per-grid in case.py to keep +# cells square; x is periodic with no gradient, so its size is physically free). +H = 1.0 # wall-normal gap (y extent) +U = 1.0 # top-wall speed (x-component) + +# Wall temperatures (bottom fixed at T0, top moving at T1) +T0 = 300.0 +T1 = 400.0 +T_ref = 0.5 * (T0 + T1) # mid-gap reference temperature + +# Soft sound speed: c = sqrt(gam*(gam-1)*cv*T). Pick c_ref at T_ref, back out cv. +c_ref = 10.0 +cv = c_ref**2 / (gam * (gam - 1.0) * T_ref) +cp = gam * cv + +# Uniform pressure giving rho = rho_ref at T_ref (mechanical equilibrium: no body +# force, so a y-varying density at uniform pressure is a valid steady state). +rho_ref = 1.0 +p0 = (gam - 1.0) * rho_ref * cv * T_ref + +# Arrhenius viscosity mu(T) = exp(C + D/T). +mu_contrast = 3.5 # target mu(T0)/mu(T1) across the gap +D = np.log(mu_contrast) / (1.0 / T0 - 1.0 / T1) +mu_ref = 0.05 # dynamic viscosity at T_ref -> Re = rho_ref*U*H/mu_ref +C = np.log(mu_ref) - D / T_ref + +# Thermal conductivity from a chosen Prandtl number Pr = mu*cp/k. +Pr = 1.0 +k_therm = mu_ref * cp / Pr + + +def mu_of_T(T): + """Arrhenius dynamic viscosity, same formula MFC evaluates (visc_model = 1).""" + return np.exp(C + D / T) + + +def T_linear(y): + """Pure-conduction (Brinkman -> 0) temperature profile across the gap.""" + return T0 + (T1 - T0) * (y / H) + + +def rho_of_T(T): + """Density that encodes temperature T at the uniform pressure p0 via the EOS.""" + return (p0 + p_inf) / ((gam - 1.0) * cv * T) + + +# Derived dimensionless diagnostics (reported, not tuned). +Re = rho_ref * U * H / mu_ref # Reynolds number at T_ref +Ma = U / c_ref # Mach number at T_ref +Br = mu_ref * U**2 / (k_therm * (T1 - T0)) # Brinkman: viscous-heating strength + +# Slowest viscous relaxation mode sets how long to run to reach steady state. +nu_min = mu_of_T(T1) / rho_of_T(T1) # smallest kinematic viscosity (hot wall) +nu_max = mu_of_T(T0) / rho_of_T(T0) # largest kinematic viscosity (cold wall) +tau_relax = H**2 / (np.pi**2 * nu_min) # first-mode decay time +t_end = 10.0 * tau_relax # run ~10 decay times: steady transient well below grid error + + +if __name__ == "__main__": + print(f"cv = {cv:.6g} cp = {cp:.6g}") + print(f"p0 = {p0:.6g}") + print(f"C = {C:.6g} D = {D:.6g}") + print(f"k_therm = {k_therm:.6g}") + print(f"mu(T0) = {mu_of_T(T0):.6g} mu(T1) = {mu_of_T(T1):.6g} contrast = {mu_of_T(T0) / mu_of_T(T1):.3f}") + print(f"Re = {Re:.4g} Ma = {Ma:.4g} Br = {Br:.4g} Pr = {Pr:.4g}") + print(f"nu_min = {nu_min:.4g} nu_max = {nu_max:.4g}") + print(f"tau_relax = {tau_relax:.4g} t_end = {t_end:.4g}") diff --git a/examples/2D_Couette_Variable_Viscosity/reference.py b/examples/2D_Couette_Variable_Viscosity/reference.py new file mode 100644 index 0000000000..dc4792b5c7 --- /dev/null +++ b/examples/2D_Couette_Variable_Viscosity/reference.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +"""Exact steady solution of variable-viscosity planar Couette flow. + +The incompressible steady state satisfies two coupled ODEs across the gap +0 <= y <= H: + + momentum: d/dy [ mu(T) du/dy ] = 0 (uniform shear stress) + energy: d/dy [ k dT/dy ] + mu(T) (du/dy)^2 = 0 (conduction + heating) + +with mu(T) = exp(C + D/T), and boundary conditions + + u(0) = 0, u(H) = U, T(0) = T0, T(H) = T1. + +Writing the momentum flux tau = mu du/dy and heat flux q = k dT/dy as states +turns this into a first-order BVP that scipy.solve_bvp resolves to ~1e-10: + + u' = tau / mu(T) + tau' = 0 + T' = q / k + q' = -tau^2 / mu(T) (= -mu (du/dy)^2, the dissipation source) + +The Brinkman-number -> 0 limit (negligible viscous heating) has the closed form +T(y) linear and u(y) = U * I(y)/I(H) with I(y) = integral_0^y dy'/mu(T(y')); we +use it as an independent check on the BVP solver. All physical constants come +from couette_config.py so this never drifts from the MFC case. +""" + +import couette_config as cfg +import numpy as np +from scipy.integrate import solve_bvp + + +def solve_couette(npts=2001): + """Solve the coupled BVP. Returns (y, u, T, tau, q_flux) on a fine mesh.""" + H, U, T0, T1, k = cfg.H, cfg.U, cfg.T0, cfg.T1, cfg.k_therm + + def rhs(y, s): + u, tau, T, q = s + mu = cfg.mu_of_T(T) + return np.vstack([tau / mu, np.zeros_like(y), q / k, -(tau**2) / mu]) + + def bc(sa, sb): + # sa = state at y=0, sb = state at y=H + return np.array([sa[0] - 0.0, sb[0] - U, sa[2] - T0, sb[2] - T1]) + + y = np.linspace(0.0, H, 51) + # Initial guess: linear u and T; stress and heat flux from those slopes. + mu_mid = cfg.mu_of_T(cfg.T_ref) + s0 = np.vstack( + [ + U * y / H, + np.full_like(y, mu_mid * U / H), + T0 + (T1 - T0) * y / H, + np.full_like(y, k * (T1 - T0) / H), + ] + ) + sol = solve_bvp(rhs, bc, y, s0, tol=1e-10, max_nodes=200000) + if not sol.success: + raise RuntimeError(f"solve_bvp failed: {sol.message}") + + yf = np.linspace(0.0, H, npts) + u, tau, T, q = sol.sol(yf) + return yf, u, T, float(tau[0]), q + + +def u_decoupled(y): + """Brinkman -> 0 closed form: linear T, u = U * I(y)/I(H).""" + yg = np.linspace(0.0, cfg.H, 20001) + integrand = 1.0 / cfg.mu_of_T(cfg.T_linear(yg)) + cum = np.concatenate([[0.0], np.cumsum(0.5 * (integrand[1:] + integrand[:-1]) * np.diff(yg))]) + return cfg.U * np.interp(y, yg, cum) / cum[-1] + + +def u_at(y): + """Exact (coupled-BVP) velocity sampled at arbitrary y.""" + yf, u, _T, _tau, _q = solve_couette() + return np.interp(y, yf, u) + + +def T_at(y): + """Exact (coupled-BVP) temperature sampled at arbitrary y.""" + yf, _u, T, _tau, _q = solve_couette() + return np.interp(y, yf, T) + + +if __name__ == "__main__": + yf, u, T, tau, q = solve_couette() + # Cross-check the coupled solver against the Br->0 closed form. With the tiny + # Brinkman number of this case the two velocity profiles agree closely; the + # gap is exactly the viscous-heating contribution. + ud = u_decoupled(yf) + rel_u_gap = np.max(np.abs(u - ud)) / cfg.U + Tlin = cfg.T_linear(yf) + print(f"shear stress tau = {tau:.6g}") + print(f"max |u_coupled - u_decoup| = {rel_u_gap * cfg.U:.3e} ({100 * rel_u_gap:.3f}% of U; ~ Brinkman effect, Br={cfg.Br:.3g})") + print(f"max viscous-heating bump = {np.max(T - Tlin):.4g} K above the linear profile") + print(f"u(H/2) curved value = {np.interp(cfg.H / 2, yf, u):.6g} (straight-line Couette would give {cfg.U / 2:.6g})") diff --git a/examples/2D_Couette_Variable_Viscosity/run_suite.sh b/examples/2D_Couette_Variable_Viscosity/run_suite.sh new file mode 100755 index 0000000000..ce82ad626a --- /dev/null +++ b/examples/2D_Couette_Variable_Viscosity/run_suite.sh @@ -0,0 +1,56 @@ +#!/usr/bin/env bash +# Run the variable-viscosity Couette grid-refinement sweep. +# +# The first grid is a normal building run (the analytic density IC compiles into +# pre_process, so the build is case-specific); once it finishes, the remaining +# grids reuse that binary via --no-build and run concurrently, each in its own +# directory (runs/n/) on a single rank so they cannot collide and never touch +# the MPI domain decomposition. Pass grid sizes as arguments; default is +# "32 64 96". +# +# ./run_suite.sh # 32 64 96 +# ./run_suite.sh 32 64 96 128 # add a finer grid (n=128 is viscous-CFL +# # limited: ~221k steps, hours on one rank) +# +# After this finishes, analyze with: python3 validate.py +set -euo pipefail + +HERE="$(cd "$(dirname "$0")" && pwd)" +ROOT="$(cd "$HERE/../.." && pwd)" +GRIDS="${*:-32 64 96}" + +pids=() +first=1 +for N in $GRIDS; do + d="$HERE/runs/n$N" + rm -rf "$d" + mkdir -p "$d" + cp "$HERE/case.py" "$HERE/couette_config.py" "$d/" + if [ "$first" -eq 1 ]; then + first=0 + echo "running n=$N -> $d (builds MFC, then the other grids launch)" + ( + cd "$ROOT" + COUETTE_N="$N" ./mfc.sh run "$d/case.py" -n 1 >"$d/run.log" 2>&1 + ) + continue + fi + echo "launching n=$N -> $d" + ( + cd "$ROOT" + COUETTE_N="$N" ./mfc.sh run "$d/case.py" -n 1 --no-build >"$d/run.log" 2>&1 + ) & + pids+=("$!") + sleep 5 # stagger so the per-run input-file generation does not overlap +done + +echo "waiting on ${#pids[@]} runs: ${pids[*]}" +rc=0 +for p in "${pids[@]}"; do + wait "$p" || rc=1 +done +if [ "$rc" -ne 0 ]; then + echo "ERROR: at least one grid failed; check runs/n*/run.log" >&2 + exit 1 +fi +echo "all grids done" diff --git a/examples/2D_Couette_Variable_Viscosity/validate.py b/examples/2D_Couette_Variable_Viscosity/validate.py new file mode 100644 index 0000000000..39a89c3f79 --- /dev/null +++ b/examples/2D_Couette_Variable_Viscosity/validate.py @@ -0,0 +1,208 @@ +#!/usr/bin/env python3 +"""Compare MFC's variable-viscosity Couette runs to the exact solution. + +This script does NOT run MFC. Run the cases first with run_suite.sh (which leaves +each grid in runs/n/), then run this to analyze them: + + ./run_suite.sh 32 64 96 + python3 validate.py # analyzes runs/n32, runs/n64, runs/n96 + +For each grid it reads back the steady velocity and temperature profiles MFC +saved, compares them to the exact coupled-BVP solution from reference.py, reports +the error, and -- across grids -- the observed spatial order of accuracy. It +writes figures/ and summary.json. + +The headline check is the velocity profile: mu(T) = exp(C + D/T) makes the hot +(thin) fluid near the moving wall shear faster than the cold (viscous) fluid near +the fixed wall, so u(y) is curved. A constant-viscosity solver would give a +straight line; the curvature is the mu(T) signal, and it must match the exact +profile to the discretization error. +""" + +import glob +import json +import os +import re +import sys + +import matplotlib +import numpy as np + +matplotlib.use("Agg") +import couette_config as cfg +import matplotlib.pyplot as plt +import reference + +HERE = os.path.dirname(os.path.abspath(__file__)) +RUNS = os.path.join(HERE, "runs") +FIG = os.path.join(HERE, "figures") +SUMMARY = os.path.join(HERE, "summary.json") + +# model_eqns = 3, num_fluids = 1, 2D conserved-field layout (0-based): +# 0: rho 1: rho*u 2: rho*v 3: E 4: alpha 5: rho*e (internal energy) +I_RHO, I_MOMX, I_RHOE = 0, 1, 5 + + +def read_inp(rundir): + inp = {} + for line in open(os.path.join(rundir, "simulation.inp")): + if "=" in line: + k, v = line.split("=", 1) + inp[k.strip().lower()] = v.strip().rstrip(",") + return inp + + +def steps_of(rundir): + rd = os.path.join(rundir, "restart_data") + steps = [] + for f in glob.glob(os.path.join(rd, "lustre_*.dat")): + match = re.search(r"lustre_(\d+)\.dat$", os.path.basename(f)) + if match: + steps.append(int(match.group(1))) + return sorted(steps) + + +def load_step(rundir, m, n, step): + a = np.fromfile(os.path.join(rundir, "restart_data", f"lustre_{step}.dat"), np.float64) + return a.reshape(a.size // (m * n), 1, n, m) # (nvars, p=1, n, m), x innermost + + +def profiles(fields): + """x-averaged wall-normal profiles u(y) and T(y) from the conserved fields.""" + rho = fields[I_RHO, 0] # (n, m) + u = fields[I_MOMX, 0] / rho + T = fields[I_RHOE, 0] / (rho * cfg.cv) # ideal gas (p_inf = 0): T = e/cv + return u.mean(axis=1), T.mean(axis=1) # average over x -> (n,) + + +def analyze_grid(N): + rundir = os.path.join(RUNS, f"n{N}") + if not os.path.isdir(os.path.join(rundir, "restart_data")): + raise SystemExit(f"missing {rundir}/restart_data -- run ./run_suite.sh {N} first") + inp = read_inp(rundir) + m, n = int(inp["m"]) + 1, int(inp["n"]) + 1 + yc = (np.arange(n) + 0.5) * cfg.H / n + steps = steps_of(rundir) + + u_last, T_last = profiles(load_step(rundir, m, n, steps[-1])) + u_prev, _ = profiles(load_step(rundir, m, n, steps[-2])) + # Steadiness: how much the velocity profile still moves between the last two + # saved snapshots, relative to the wall speed. + unsteadiness = float(np.max(np.abs(u_last - u_prev)) / cfg.U) + + u_ex = reference.u_at(yc) + T_ex = reference.T_at(yc) + err_u = float(np.sqrt(np.mean((u_last - u_ex) ** 2)) / cfg.U) + err_u_inf = float(np.max(np.abs(u_last - u_ex)) / cfg.U) + err_T = float(np.sqrt(np.mean((T_last - T_ex) ** 2)) / (cfg.T1 - cfg.T0)) + return { + "N": N, + "ny": n, + "dy": cfg.H / n, + "nsteps": steps[-1], + "unsteadiness": unsteadiness, + "err_u_L2": err_u, + "err_u_Linf": err_u_inf, + "err_T_L2": err_T, + "yc": yc, + "u": u_last, + "T": T_last, + "u_ex": u_ex, + "T_ex": T_ex, + } + + +def main(grids): + os.makedirs(FIG, exist_ok=True) + results = [analyze_grid(N) for N in grids] + results.sort(key=lambda r: r["ny"]) + + nys = np.array([r["ny"] for r in results]) + errs_u = np.array([r["err_u_L2"] for r in results]) + errs_T = np.array([r["err_T_L2"] for r in results]) + # error ~ dy^order; dy = H/Ny, so fit log(err) vs log(dy) -> positive order + dys = cfg.H / nys + order_u = float(np.polyfit(np.log(dys), np.log(errs_u), 1)[0]) + + # --- console report --- + print(f"\nVariable-viscosity Couette validation (Re={cfg.Re:.0f}, Ma={cfg.Ma:.2f}, Br={cfg.Br:.3f}, mu-contrast={cfg.mu_of_T(cfg.T0) / cfg.mu_of_T(cfg.T1):.2f})") + print(f"{'Ny':>5} {'dy':>10} {'unsteady':>11} {'L2(u)/U':>11} {'Linf(u)/U':>11} {'L2(T)/dT':>11}") + for r in results: + print(f"{r['ny']:>5} {r['dy']:>10.3e} {r['unsteadiness']:>11.2e} {r['err_u_L2']:>11.3e} {r['err_u_Linf']:>11.3e} {r['err_T_L2']:>11.3e}") + print(f"observed spatial order (L2 velocity) = {order_u:.2f}") + + # --- figure 1: profile overlay on the finest grid --- + rf = results[-1] + fig, (axu, axt) = plt.subplots(1, 2, figsize=(11, 4.6)) + axu.plot(rf["u_ex"], rf["yc"], "-", color="k", lw=2, label="exact (coupled BVP)") + axu.plot(rf["u"], rf["yc"], "o", ms=4, color="C0", label=f"MFC (Ny={rf['ny']})") + axu.plot(cfg.U * rf["yc"] / cfg.H, rf["yc"], "--", color="C3", lw=1, label=r"constant-$\mu$ (straight)") + axu.set_xlabel("u / U"), axu.set_ylabel("y / H") + axu.set_title("velocity: curvature is the mu(T) signal"), axu.legend(fontsize=8) + axt.plot(rf["T_ex"], rf["yc"], "-", color="k", lw=2, label="exact") + axt.plot(rf["T"], rf["yc"], "s", ms=4, color="C1", label="MFC") + axt.set_xlabel("T [K]"), axt.set_ylabel("y / H") + axt.set_title("temperature profile"), axt.legend(fontsize=8) + fig.tight_layout() + fig.savefig(os.path.join(FIG, "couette_profiles.png"), dpi=130) + + # --- figure 2: convergence --- + fig2, ax = plt.subplots(figsize=(6, 5)) + ax.loglog(nys, errs_u, "o-", label=f"L2 velocity (order {order_u:.2f})") + ax.loglog(nys, errs_T, "s-", label="L2 temperature") + ref2 = errs_u[0] * (nys[0] / nys) ** 2 + ax.loglog(nys, ref2, "k--", lw=1, label="2nd-order reference") + ax.set_xlabel("Ny (wall-normal cells)"), ax.set_ylabel("relative L2 error") + ax.set_title("grid convergence, mu(T) Couette"), ax.legend() + ax.grid(True, which="both", alpha=0.3) + fig2.tight_layout() + fig2.savefig(os.path.join(FIG, "couette_convergence.png"), dpi=130) + + # --- figure 3: spatial relative error --- + # Left: signed relative velocity error along the gap, one curve per grid -- + # shows where the error lives and that it shrinks under refinement. Right: the + # same error scaled by (Ny/Ny_coarsest)^2; the curves collapse onto one shape, + # so the 2nd-order convergence holds pointwise, not just in the L2 norm. + ny0 = nys[0] + fig3, (axe, axs) = plt.subplots(1, 2, figsize=(11, 4.6)) + for r in results: + rel = 1e4 * (r["u"] - r["u_ex"]) / cfg.U # in units of 1e-4 for readable ticks + axe.plot(rel, r["yc"], "-o", ms=3, label=f"Ny={r['ny']}") + axs.plot(rel * (r["ny"] / ny0) ** 2, r["yc"], "-o", ms=3, label=f"Ny={r['ny']}") + for a in (axe, axs): + a.axvline(0.0, color="k", lw=0.6) + a.set_ylabel("y / H") + a.legend(fontsize=8) + a.xaxis.set_major_locator(plt.MaxNLocator(5)) + axe.set_xlabel(r"$(u_{\rm MFC} - u_{\rm exact})\,/\,U\quad[\times 10^{-4}]$") + axe.set_title("relative velocity error") + axs.set_xlabel(r"$(u_{\rm MFC} - u_{\rm exact})/U \times (N_y/N_{y,0})^2\quad[\times 10^{-4}]$") + axs.set_title("grid-scaled error (collapse = 2nd order in space)") + fig3.tight_layout() + fig3.savefig(os.path.join(FIG, "couette_error.png"), dpi=130) + + # --- summary.json --- + summary = { + "config": { + "Re": cfg.Re, + "Ma": cfg.Ma, + "Br": cfg.Br, + "Pr": cfg.Pr, + "mu_contrast": cfg.mu_of_T(cfg.T0) / cfg.mu_of_T(cfg.T1), + "T0": cfg.T0, + "T1": cfg.T1, + "U": cfg.U, + "u_mid_exact": float(reference.u_at(np.array([cfg.H / 2]))[0]), + "u_mid_constant_mu": cfg.U / 2, + }, + "observed_order_u_L2": order_u, + "grids": [{k: r[k] for k in ("N", "ny", "dy", "nsteps", "unsteadiness", "err_u_L2", "err_u_Linf", "err_T_L2")} for r in results], + } + json.dump(summary, open(SUMMARY, "w"), indent=2) + print(f"\nwrote {SUMMARY}") + print(f"wrote {FIG}/couette_profiles.png, couette_convergence.png, couette_error.png") + + +if __name__ == "__main__": + grids = [int(a) for a in sys.argv[1:]] or [32, 64, 96] + main(grids) diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index 4ceb53b887..8643a89b63 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -3327,6 +3327,7 @@ def foreach_example(): "1D_conduction_convergence", "2D_axisym_conduction_convergence", "3D_cyl_azimuthal_conduction_convergence", + "2D_Couette_Variable_Viscosity", # analytic IC; mu(T) is covered by the "Viscous -> visc_model=1" case "3D_IGR_TaylorGreenVortex_nvidia", "2D_backward_facing_step", "2D_forward_facing_step", From 804ed64d8bcae4c32581edfa002e2f7b6596e017 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Fri, 23 Oct 2026 16:53:49 -0400 Subject: [PATCH 3/4] Pass mixture temperature into f_compute_arrhenius_shear_re nvfortran 25.5 crashes (fort2 signal 11) inlining a seq GPU routine that itself calls the inlined f_mixture_temperature into the m_data_output and m_time_steppers kernels under --case-optimization --gpu acc/mp. Callers now compute T with f_mixture_temperature and pass it in, so the routine has no nested inlined call. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VtsGVvCfsDJz5dmtSk9dKt --- src/simulation/m_arrhenius_viscosity.fpp | 17 +++++++++-------- src/simulation/m_data_output.fpp | 7 +++++-- src/simulation/m_riemann_solver_hll.fpp | 6 ++++-- src/simulation/m_riemann_solver_hllc.fpp | 12 ++++++++---- src/simulation/m_time_steppers.fpp | 4 +++- 5 files changed, 29 insertions(+), 17 deletions(-) diff --git a/src/simulation/m_arrhenius_viscosity.fpp b/src/simulation/m_arrhenius_viscosity.fpp index 1c56b232ce..35f5185f6b 100644 --- a/src/simulation/m_arrhenius_viscosity.fpp +++ b/src/simulation/m_arrhenius_viscosity.fpp @@ -10,7 +10,6 @@ module m_arrhenius_viscosity use m_derived_types !< Definitions of the derived types use m_global_parameters !< sgm_eps, visc_models, visc_cs, visc_ds, fluid_inv_re use m_constants, only: visc_exp_arg_cap - use m_eos, only: f_mixture_temperature implicit none @@ -18,23 +17,25 @@ module m_arrhenius_viscosity contains - !> Mixture shear Reynolds number 1/sum(alpha_i*mu_i) of one state at its thermal-equilibrium mixture temperature, the same - !! closure heat conduction uses. Arrhenius fluids take mu_i = exp(C_i + D_i/T); the others keep their constant 1/Re_i(1). - function f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) result(Re_shear) + !> Mixture shear Reynolds number 1/sum(alpha_i*mu_i) of one state at its mixture temperature T_mix, which callers take from + !! f_mixture_temperature (the closure heat conduction uses). Arrhenius fluids take mu_i = exp(C_i + D_i/T); the others keep + !! their constant 1/Re_i(1). T_mix is an argument rather than computed here because nvfortran 25.5 crashes (fort2 signal 11) + !! inlining the nested f_mixture_temperature call into the m_data_output and m_time_steppers kernels. + function f_compute_arrhenius_shear_re(alpha, T_mix) result(Re_shear) $:GPU_ROUTINE(function_name='f_compute_arrhenius_shear_re', parallelism='[seq]', cray_inline=True) #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha, alpha_rho + real(wp), dimension(3), intent(in) :: alpha #:else - real(wp), dimension(num_fluids), intent(in) :: alpha, alpha_rho + real(wp), dimension(num_fluids), intent(in) :: alpha #:endif - real(wp), intent(in) :: pres, gamma, pi_inf + real(wp), intent(in) :: T_mix real(wp) :: Re_shear real(wp) :: T, inv_re integer :: fl - T = max(f_mixture_temperature(alpha_rho, pres, gamma, pi_inf), sgm_eps) + T = max(T_mix, sgm_eps) inv_re = 0._wp do fl = 1, num_fluids diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 01f5ae1f0d..5d453c243f 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -17,6 +17,7 @@ module m_data_output use m_helper_basic use m_sim_helpers use m_arrhenius_viscosity + use m_eos, only: f_mixture_temperature use m_delay_file_access use m_ibm use m_boundary_common @@ -238,7 +239,8 @@ contains Re(1) = 1._wp/max(Re(1), sgm_eps) end if - if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) + if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, f_mixture_temperature(alpha_rho, & + & pres, gamma, pi_inf)) call s_compute_stability_from_dt(vel, c, rho, Re, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) @@ -401,7 +403,8 @@ contains Re(1) = 1._wp/max(Re(1), sgm_eps) end if - if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) + if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, f_mixture_temperature(alpha_rho, pres, & + & gamma, pi_inf)) call s_compute_stability_from_dt(vel, c, rho, Re, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 0fe0e675df..fa03075710 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -202,8 +202,10 @@ contains call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) call s_compute_interface_reynolds(alpha_R, Re_R, Re_size_loc1, Re_size_loc2) if (any_arrhenius_visc) then - Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, alpha_rho_L, pres_L, gamma_L, pi_inf_L) - Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, alpha_rho_R, pres_R, gamma_R, pi_inf_R) + Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, f_mixture_temperature(alpha_rho_L, pres_L, & + & gamma_L, pi_inf_L)) + Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, f_mixture_temperature(alpha_rho_R, pres_R, & + & gamma_R, pi_inf_R)) end if end if diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index f6a3d814b0..be864c0f8c 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -258,8 +258,10 @@ contains call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) call s_compute_interface_reynolds(alpha_R, Re_R, Re_size_loc1, Re_size_loc2) if (any_arrhenius_visc) then - Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, alpha_rho_L, pres_L, gamma_L, pi_inf_L) - Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, alpha_rho_R, pres_R, gamma_R, pi_inf_R) + Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, f_mixture_temperature(alpha_rho_L, & + & pres_L, gamma_L, pi_inf_L)) + Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, f_mixture_temperature(alpha_rho_R, & + & pres_R, gamma_R, pi_inf_R)) end if end if @@ -998,8 +1000,10 @@ contains call s_compute_interface_reynolds(alpha_L, Re_L, Re_size_loc1, Re_size_loc2) call s_compute_interface_reynolds(alpha_R, Re_R, Re_size_loc1, Re_size_loc2) if (any_arrhenius_visc) then - Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, alpha_rho_L, pres_L, gamma_L, pi_inf_L) - Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, alpha_rho_R, pres_R, gamma_R, pi_inf_R) + Re_L(1) = f_compute_arrhenius_shear_re(alpha_L, f_mixture_temperature(alpha_rho_L, & + & pres_L, gamma_L, pi_inf_L)) + Re_R(1) = f_compute_arrhenius_shear_re(alpha_R, f_mixture_temperature(alpha_rho_R, & + & pres_R, gamma_R, pi_inf_R)) end if end if diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index c2b4983242..ec9b20810c 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -26,6 +26,7 @@ module m_time_steppers use m_helper use m_sim_helpers use m_arrhenius_viscosity + use m_eos, only: f_mixture_temperature use m_fftw use m_nvtx use m_thermochem, only: num_species @@ -734,7 +735,8 @@ contains Re(1) = 1._wp/max(Re(1), sgm_eps) end if - if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, alpha_rho, pres, gamma, pi_inf) + if (any_arrhenius_visc) Re(1) = f_compute_arrhenius_shear_re(alpha, f_mixture_temperature(alpha_rho, & + & pres, gamma, pi_inf)) call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, alpha, alpha_rho, j, k, l) From 0e1527e071db7569ed96f91b3861e2945ed6802d Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Fri, 23 Oct 2026 17:42:44 -0400 Subject: [PATCH 4/4] Include case.fpp in m_arrhenius_viscosity Without it MFC_CASE_OPTIMIZATION keeps its CMake default (False) in this module, so an amdflang --case-optimization build declared the alpha dummy dimension(3) while HLL/HLLC pass dimension(num_fluids) actuals, which amdflang rejects ("Actual argument array has fewer elements (2) than dummy argument 'alpha=' array (3)", Frontier AMD case-opt gpu-omp). Other compilers' preprocessed output is unchanged apart from comments. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VtsGVvCfsDJz5dmtSk9dKt --- src/simulation/m_arrhenius_viscosity.fpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/simulation/m_arrhenius_viscosity.fpp b/src/simulation/m_arrhenius_viscosity.fpp index 35f5185f6b..d7c767bdb6 100644 --- a/src/simulation/m_arrhenius_viscosity.fpp +++ b/src/simulation/m_arrhenius_viscosity.fpp @@ -2,6 +2,7 @@ !! @file m_arrhenius_viscosity.fpp !! @brief Temperature-dependent (Arrhenius) shear viscosity mu = exp(C + D/T) and the mixture shear Reynolds number it gives. +#:include 'case.fpp' #:include 'macros.fpp' !> @brief Temperature-dependent (Arrhenius) shear viscosity mu = exp(C + D/T) and the mixture shear Reynolds number it gives.