diff --git a/examples/1D_mg_acoustic/case.py b/examples/1D_mg_acoustic/case.py new file mode 100644 index 000000000..c14071434 --- /dev/null +++ b/examples/1D_mg_acoustic/case.py @@ -0,0 +1,70 @@ +""" +Right-moving acoustic pulse in a single Mie-Gruneisen fluid at its reference state. +The pulse is a simple wave (drho, dp = c^2 drho, du = c drho/rho0) so only the right-going +characteristic carries it; the harness measures its speed against the general analytic c. +""" + +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="1D Mie-Gruneisen acoustic pulse") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=200) +parser.add_argument("--cfl", type=float, default=0.4) +args = parser.parse_args() + +rho0, p0, c0, s, gruneisen = 1.0, 1.0, 1.0, 1.5, 0.4 +c = math.sqrt(c0**2 + (1.0 + gruneisen) * p0 / rho0) # the frozen speed at the reference state +amp, x0, width = 1.0e-4, 0.3, 0.05 +N, L, T_end = args.N, 1.0, 0.4 +dt = args.cfl * (L / N) / c +Nt = math.ceil(T_end / dt) +dt = T_end / Nt +pulse = f"{amp}*exp(-((x - {x0})/{width})**2)" + +print( + json.dumps( + { + "run_time_info": "F", + "x_domain%beg": 0.0, + "x_domain%end": L, + "m": N - 1, + "n": 0, + "p": 0, + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "num_patches": 1, + "model_eqns": 2, + "num_fluids": 1, + "time_stepper": 3, + "recon_type": "weno", + "weno_order": 5, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%length_x": L, + "patch_icpp(1)%alpha_rho(1)": f"{rho0} + {pulse}", + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%vel(1)": f"{c}/{rho0}*{pulse}", + "patch_icpp(1)%pres": f"{p0} + {c}**2*{pulse}", + "fluid_pp(1)%eos": "mie_gruneisen", + "fluid_pp(1)%mg_rho0": rho0, + "fluid_pp(1)%mg_c0": c0, + "fluid_pp(1)%mg_s": s, + "fluid_pp(1)%mg_gruneisen": gruneisen, + } + ) +) diff --git a/examples/1D_mg_impact/case.py b/examples/1D_mg_impact/case.py new file mode 100644 index 000000000..619a966a7 --- /dev/null +++ b/examples/1D_mg_impact/case.py @@ -0,0 +1,71 @@ +""" +Symmetric impact of two Mie-Gruneisen slabs approaching at relative speed U. +Each slab is brought to rest by a shock with particle-velocity jump U/2, so the shock state +lies exactly on the Hugoniot u_s = c0 + s u_p; the harness checks the shock speed and the +plateau density against that relation. +""" + +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="1D Mie-Gruneisen symmetric impact") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=800) +parser.add_argument("--U", type=float, default=1.0, help="closing speed of the two slabs") +parser.add_argument("--cfl", type=float, default=0.4) +args = parser.parse_args() + +rho0, p0, c0, s, gruneisen = 1.0, 1.0e-3, 1.0, 1.5, 0.4 +N, L, T_end = args.N, 1.0, 0.2 +dt = args.cfl * (L / N) / (c0 + (s + 1.0) * args.U) +Nt = math.ceil(T_end / dt) +dt = T_end / Nt + +case = { + "run_time_info": "F", + "x_domain%beg": 0.0, + "x_domain%end": L, + "m": N - 1, + "n": 0, + "p": 0, + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "num_patches": 2, + "model_eqns": 2, + "num_fluids": 1, + "time_stepper": 3, + "recon_type": "weno", + "weno_order": 5, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "fluid_pp(1)%eos": "mie_gruneisen", + "fluid_pp(1)%mg_rho0": rho0, + "fluid_pp(1)%mg_c0": c0, + "fluid_pp(1)%mg_s": s, + "fluid_pp(1)%mg_gruneisen": gruneisen, +} +for pid, (x_c, vel) in enumerate([(0.25, 0.5 * args.U), (0.75, -0.5 * args.U)], start=1): + case.update( + { + f"patch_icpp({pid})%geometry": 1, + f"patch_icpp({pid})%x_centroid": x_c, + f"patch_icpp({pid})%length_x": 0.5 * L, + f"patch_icpp({pid})%alpha_rho(1)": rho0, + f"patch_icpp({pid})%alpha(1)": 1.0, + f"patch_icpp({pid})%vel(1)": vel, + f"patch_icpp({pid})%pres": p0, + } + ) +print(json.dumps(case)) diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 5f9088961..f317ef9bf 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -118,6 +118,7 @@ module m_constants !! cannot be auto-generated, so these are hand-written. integer, parameter :: eos_stiffened_gas = 1 integer, parameter :: eos_ideal_gas = 2 + integer, parameter :: eos_mie_gruneisen = 3 integer, parameter :: num_synth_shells_max = 50 !< Max energy shells for synthetic turbulence integer, parameter :: num_turb_sources_max = 10 !< Max Gaussian forcing zones for synthetic turbulence diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index d80c885b4..bca9ca5ee 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -394,6 +394,10 @@ module m_derived_types real(wp) :: qvp !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004)) real(wp) :: G integer :: eos !< Equation of state selector (eos_* in m_constants) + real(wp) :: mg_rho0 !< Mie-Gruneisen reference density + real(wp) :: mg_c0 !< Mie-Gruneisen bulk sound speed at mg_rho0 + real(wp) :: mg_s !< Mie-Gruneisen linear Hugoniot slope, u_s = c0 + s u_p + real(wp) :: mg_gruneisen !< Gruneisen coefficient Gamma_G (not the shear modulus G) logical :: non_newtonian !< Enable Herschel-Bulkley non-Newtonian viscosity real(wp) :: K !< HB consistency index real(wp) :: nn !< HB flow behavior index diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 3c230ff66..7c46f1458 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -54,6 +54,11 @@ module m_global_parameters_common !! written as p + B = const*rho**n. real(wp), allocatable, dimension(:) :: gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps $:GPU_DECLARE(create='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps]') + !> Per-fluid EOS selector and Mie-Gruneisen reference curve, resolved once at init like the arrays above + integer, allocatable, dimension(:) :: eoss + real(wp), allocatable, dimension(:) :: mg_rho0s, mg_c0s, mg_ss, mg_gruneisens + logical :: any_state_dependent_eos !< True when some fluid's coefficients vary with density; set at init + $:GPU_DECLARE(create='[eoss, mg_rho0s, mg_c0s, mg_ss, mg_gruneisens, any_state_dependent_eos]') !> @} !> @name Fluids participating in shear and bulk viscosity diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index f8bbd640f..75e548c54 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -28,8 +28,8 @@ module m_variables_conversion & s_compute_species_fraction, s_compute_mixture_coefficients, s_compute_energy, s_compute_speed_of_sound, f_bulk_modulus, & & f_pressure, f_phase_internal_energy, f_isentrope_exponent, f_isentrope_pressure, f_sg_thermal, f_pressure_on_isentrope, & & s_compute_mixture_coefficients_dt, s_compute_speed_of_sound_avg, s_compute_fast_magnetosonic_speed, f_elastic_energy, & - & f_hypoelastic_energy, f_relativistic_enthalpy, s_finalize_variables_conversion_module, gammas, isentrope_n, pi_infs, & - & isentrope_B, cvs, qvs, qvps + & f_hypoelastic_energy, f_relativistic_enthalpy, s_eos_coefficients, s_finalize_variables_conversion_module, gammas, & + & isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps real(wp), allocatable, dimension(:) :: Gs_vc integer, allocatable, dimension(:) :: bubrs_vc @@ -267,6 +267,8 @@ contains $:GPU_UPDATE(device='[enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc]') @:ALLOCATE(gammas (1:num_fluids)) + @:ALLOCATE(eoss (1:num_fluids), mg_rho0s (1:num_fluids), mg_c0s (1:num_fluids), mg_ss (1:num_fluids), & + & mg_gruneisens (1:num_fluids)) @:ALLOCATE(isentrope_n (1:num_fluids)) @:ALLOCATE(pi_infs(1:num_fluids)) @:ALLOCATE(isentrope_B(1:num_fluids)) @@ -275,6 +277,7 @@ contains @:ALLOCATE(qvps (1:num_fluids)) @:ALLOCATE(Gs_vc (1:num_fluids)) + any_state_dependent_eos = .false. do i = 1, num_fluids gammas(i) = fluid_pp(i)%gamma isentrope_n(i) = f_isentrope_exponent(gammas(i)) @@ -292,8 +295,15 @@ contains cvs(i) = fluid_pp(i)%cv qvs(i) = fluid_pp(i)%qv qvps(i) = fluid_pp(i)%qvp + eoss(i) = fluid_pp(i)%eos + mg_rho0s(i) = fluid_pp(i)%mg_rho0 + mg_c0s(i) = fluid_pp(i)%mg_c0 + mg_ss(i) = fluid_pp(i)%mg_s + mg_gruneisens(i) = fluid_pp(i)%mg_gruneisen + if (fluid_pp(i)%eos == eos_mie_gruneisen) any_state_dependent_eos = .true. end do - $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc]') + $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, mg_rho0s, mg_c0s, mg_ss, & + & mg_gruneisens, any_state_dependent_eos]') @:ALLOCATE(Res_vc(1:2, 1:max(1, Re_size_max))) Res_vc = dflt_real @@ -1173,7 +1183,7 @@ contains if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf) - @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc) + @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, mg_rho0s, mg_c0s, mg_ss, mg_gruneisens) if (allocated(bubrs_vc)) then @:DEALLOCATE(bubrs_vc) end if @@ -1196,6 +1206,7 @@ contains real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K #:endif real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K + real(wp) :: gamma_i, pi_inf_i, dpi_i integer :: i !< Loop iterator over fluids ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients @@ -1217,8 +1228,14 @@ contains $:GPU_LOOP(parallelism='[seq]') do i = 1, num_fluids rho_K = rho_K + alpha_rho_K(i) - gamma_K = gamma_K + alpha_K(i)*gammas(i) - pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) + if (any_state_dependent_eos) then + call s_eos_coefficients(alpha_rho_K(i)/max(alpha_K(i), sgm_eps), i, gamma_i, pi_inf_i, dpi_i) + else + gamma_i = gammas(i) + pi_inf_i = pi_infs(i) + end if + gamma_K = gamma_K + alpha_K(i)*gamma_i + pi_inf_K = pi_inf_K + alpha_K(i)*pi_inf_i qv_K = qv_K + alpha_rho_K(i)*qvs(i) end do end if @@ -1284,6 +1301,49 @@ contains end subroutine s_compute_energy + !> Coefficients of fluid i at density rho in the form rho e = Gamma p + Pi that every operator here consumes, with dPi/drho. Any + !! Mie-Gruneisen EOS p = p_ref + rho Gamma_G (e - e_ref) is this form with Gamma = 1/Gamma_G and Pi = rho e_ref - p_ref/Gamma_G; + !! a new family adds one case supplying its reference curve. Gamma_G is constant, so dGamma/drho = 0. Stiffened and ideal gas + !! keep the constants resolved at init, bit for bit. + subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: rho + integer, intent(in) :: i + real(wp), intent(out) :: gamma, pi_inf, dpi + real(wp) :: mu, d, p_ref, e_ref, dp_dmu, de_dmu, G0 + + select case (eoss(i)) + case (eos_mie_gruneisen) + ! Linear-Hugoniot reference curve, u_s = c0 + s u_p: p_H = rho0 c0^2 mu (1 + mu)/(1 - (s - 1) mu)^2 on + ! compression, extended linearly on release, with the Hugoniot energy e_H = p_H mu/(2 rho0 (1 + mu)). + ! Pole at mu = 1/(s - 1); the validator warns when the initial state is near it. + mu = rho/mg_rho0s(i) - 1._wp + if (mu >= 0._wp) then + d = 1._wp - (mg_ss(i) - 1._wp)*mu + p_ref = mg_rho0s(i)*mg_c0s(i)**2*mu*(1._wp + mu)/(d*d) + dp_dmu = mg_rho0s(i)*mg_c0s(i)**2*((1._wp + 2._wp*mu)*d + 2._wp*(mg_ss(i) - 1._wp)*mu*(1._wp + mu))/(d*d*d) + else + p_ref = mg_rho0s(i)*mg_c0s(i)**2*mu + dp_dmu = mg_rho0s(i)*mg_c0s(i)**2 + end if + e_ref = p_ref*mu/(2._wp*mg_rho0s(i)*(1._wp + mu)) + de_dmu = (dp_dmu*mu*(1._wp + mu) + p_ref)/(2._wp*mg_rho0s(i)*(1._wp + mu)**2) + G0 = mg_gruneisens(i) + case default + gamma = gammas(i) + pi_inf = pi_infs(i) + dpi = 0._wp + return + end select + + gamma = 1._wp/G0 + pi_inf = rho*e_ref - p_ref/G0 + dpi = e_ref + (rho*de_dmu - dp_dmu/G0)/mg_rho0s(i) ! d/drho = (1/rho0) d/dmu + + end subroutine s_eos_coefficients + !> Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n. function f_isentrope_exponent(gamma) result(n) @@ -1420,7 +1480,7 @@ contains !> Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H, |u|^2 and qv all cancel out of c^2 !! = ((Gamma + 1)p + Pi)/(Gamma rho). Averaged states, whose enthalpy is a free input, use the _avg variant. - subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c) + subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) $:GPU_ROUTINE(parallelism='[seq]') @@ -1431,8 +1491,14 @@ contains real(wp), dimension(num_fluids), intent(in) :: adv #:endif real(wp), intent(out) :: c - real(wp) :: alf !< Subgrid void fraction; dilute by construction - integer :: q + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in), optional :: alpha_rho + #:else + real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho + #:endif + real(wp) :: alf !< Subgrid void fraction; dilute by construction + real(wp) :: rho_q, gamma_q, pi_inf_q, dpi_q + integer :: q if (chemistry) then ! Reacting mixture sound speed c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho) @@ -1441,7 +1507,16 @@ contains else ! Every case below is a bulk modulus over a density. The equation of state enters ! only through f_bulk_modulus; the cases differ in how the phases are mixed. - if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean + if (any_state_dependent_eos .and. present(alpha_rho)) then ! frozen mixing: each phase's modulus at its own density + c = 0._wp + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + rho_q = alpha_rho(q)/max(adv(q), sgm_eps) + call s_eos_coefficients(rho_q, q, gamma_q, pi_inf_q, dpi_q) + c = c + adv(q)*(f_bulk_modulus(pres, gamma_q, pi_inf_q) - rho_q*dpi_q/gamma_q) + end do + c = c/rho + else if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean c = 1._wp/(rho*(adv(1)/f_bulk_modulus(pres, gammas(1), pi_infs(1)) + adv(2)/f_bulk_modulus(pres, gammas(2), & & pi_infs(2)))) else if (model_eqns == model_eqns_6eq) then ! volume-weighted arithmetic mean @@ -1475,7 +1550,7 @@ contains !> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its !! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from !! s_compute_speed_of_sound; keep the condition below in step with the branch list there. - subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c) + subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho) $:GPU_ROUTINE(parallelism='[seq]') @@ -1486,17 +1561,23 @@ contains real(wp), dimension(num_fluids), intent(in) :: adv #:endif real(wp), intent(out) :: c + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in), optional :: alpha_rho + #:else + real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho + #:endif if (chemistry) then ! Reacting mixture sound speed if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H)) else - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) end if else if (relativity) then ! Relativistic sound speed c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H) - else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler)) then - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c) + else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) & + & .or. any_state_dependent_eos) then + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) else ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives c = (H - 5.e-1*vel_sum - qv/rho)/gamma diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 1542f161a..67f0bc3f6 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -1281,7 +1281,7 @@ contains call s_compute_mixture_coefficients(alpha_rho, adv, rho, gamma, pi_inf, qv) - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho) Ma = maxvel/c if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 34effbcb9..725c92490 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -203,6 +203,10 @@ contains ! Fluids physical parameters (post-specific; G = dflt_real differs from pre/sim) do i = 1, num_fluids_max fluid_pp(i)%eos = eos_stiffened_gas + fluid_pp(i)%mg_rho0 = dflt_real + fluid_pp(i)%mg_c0 = dflt_real + fluid_pp(i)%mg_s = dflt_real + fluid_pp(i)%mg_gruneisen = dflt_real fluid_pp(i)%gamma = dflt_real fluid_pp(i)%pi_inf = dflt_real fluid_pp(i)%cv = 0._wp diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 80aeb4573..e3d89de06 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -67,7 +67,7 @@ contains integer :: t_step real(wp), dimension(nb) :: nRtmp real(wp) :: nbub - real(wp) :: gamma, lit_gamma, pi_inf, qv + real(wp) :: gamma, pi_inf, qv real(wp) :: rho real(wp) :: pres, T real(wp) :: rhoYks(1:num_species) @@ -165,7 +165,6 @@ contains end if gamma = gammas(1) - lit_gamma = isentrope_n(1) pi_inf = pi_infs(1) qv = qvs(1) @@ -199,8 +198,6 @@ contains call s_convert_to_mixture_variables(q_cons_vf, j, 0, 0, rho, gamma, pi_inf, qv) - lit_gamma = f_isentrope_exponent(gamma) - if ((i >= eqn_idx%species%beg) .and. (i <= eqn_idx%species%end)) then write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)/rho else if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) & diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index d4a5c559b..2d75411c7 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -387,6 +387,10 @@ contains ! Fluids physical parameters do i = 1, num_fluids_max fluid_pp(i)%eos = eos_stiffened_gas + fluid_pp(i)%mg_rho0 = dflt_real + fluid_pp(i)%mg_c0 = dflt_real + fluid_pp(i)%mg_s = dflt_real + fluid_pp(i)%mg_gruneisen = dflt_real fluid_pp(i)%gamma = dflt_real fluid_pp(i)%pi_inf = dflt_real fluid_pp(i)%cv = 0._wp diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 6c923b334..3467de31b 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -163,11 +163,11 @@ contains real(wp) :: rho !< Cell-avg. density #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha !< Cell-avg. volume fraction - real(wp), dimension(3) :: vel !< Cell-avg. velocity + real(wp), dimension(3) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density + real(wp), dimension(3) :: vel !< Cell-avg. velocity #:else - real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction - real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity + real(wp), dimension(num_fluids) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density + real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity #:endif real(wp) :: vel_sum !< Cell-avg. velocity sum real(wp) :: pres !< Cell-avg. pressure @@ -189,15 +189,15 @@ contains ccfl_max_loc = 0._wp Rc_min_loc = huge(1.0_wp) ! Computing Stability Criteria at Current Time-step - $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, icfl, vcfl, & - & Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', & + $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, vel, alpha, alpha_rho, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, & + & icfl, vcfl, Rc, ccfl, fl]', reduction='[[icfl_max_loc, vcfl_max_loc, ccfl_max_loc], [Rc_min_loc]]', & & reductionOp='[max, min]') do l = 0, p do k = 0, n do j = 0, m - call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, vel, vel_sum, qv, j, k, l) + call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, k, l) - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho) if (any_non_newtonian) then Re(1) = 0._wp @@ -303,7 +303,6 @@ contains logical :: file_exist !< Logical used to check existence of current time-step directory character(LEN=15) :: FMT integer :: i, j, k, l, r - real(wp) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all' write (t_step_dir, '(a,i0,a,i0)') trim(case_dir) // '/p_all/p', proc_rank, '/', t_step @@ -377,11 +376,6 @@ contains call s_write_serial_ib_data(t_step) end if - gamma = gammas(1) - lit_gamma = isentrope_n(1) - pi_inf = pi_infs(1) - qv = qvs(1) - if (precision == precision_single) then FMT = "(2F30.3)" else diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index e1a258766..618196f07 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -444,6 +444,10 @@ contains ! Fluids physical parameters (sim-specific; Re(:) and G=0._wp differ from post) do i = 1, num_fluids_max fluid_pp(i)%eos = eos_stiffened_gas + fluid_pp(i)%mg_rho0 = dflt_real + fluid_pp(i)%mg_c0 = dflt_real + fluid_pp(i)%mg_s = dflt_real + fluid_pp(i)%mg_gruneisen = dflt_real fluid_pp(i)%gamma = dflt_real fluid_pp(i)%pi_inf = dflt_real fluid_pp(i)%cv = 0._wp diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index a6413217e..6267e074c 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -333,13 +333,13 @@ contains end if end if - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R) if (wave_speeds == wave_speeds_pressure) then call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, & - & H_avg, c_sum_Yi_Phi, alpha_R, c_avg) + & H_avg, c_sum_Yi_Phi, alpha_R, c_avg, alpha_rho_R) end if if (mhd) then diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index 6d4710a10..671e81f4e 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -276,15 +276,15 @@ contains & qv_R, rho_avg, vel_avg_rms, H_avg, gamma_avg, qv_avg) end if - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R) ! Only the pressure-based wave-speed estimate reads the averaged state, and building it ! costs eight square roots per face under the Roe average. if (wave_speeds == wave_speeds_pressure) then call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, & - & H_avg, c_sum_Yi_Phi, alpha_R, c_avg) + & H_avg, c_sum_Yi_Phi, alpha_R, c_avg, alpha_rho_R) end if if (viscous) then @@ -630,9 +630,9 @@ contains end do end if - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R) ! Only the pressure-based wave-speed estimate reads the averaged state, and building it ! costs eight square roots per face under the Roe average. @@ -640,7 +640,7 @@ contains ! Zero, not c_sum_Yi_Phi: this loop never forms the chemistry average, and ! chemistry with bubbles_euler/qbmm is prohibited, so the branch is unreachable. call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, & - & H_avg, 0._wp, alpha_R, c_avg) + & H_avg, 0._wp, alpha_R, c_avg, alpha_rho_R) end if if (viscous) then @@ -1058,15 +1058,16 @@ contains end if end if - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R) ! Only the pressure-based wave-speed estimate reads the averaged state, and building it ! costs eight square roots per face under the Roe average. if (wave_speeds == wave_speeds_pressure) then call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, & - & vel_avg_rms, H_avg, c_sum_Yi_Phi, alpha_R, c_avg) + & vel_avg_rms, H_avg, c_sum_Yi_Phi, alpha_R, c_avg, & + & alpha_rho_R) end if if (viscous) then diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp index 627466678..37822dc08 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -211,9 +211,9 @@ contains call s_compute_energy(pres_R, alpha_rho_R, alpha_R, vel_R_rms, E_R) end if - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L) + call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, alpha_L, c_L, alpha_rho_L) - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R) + call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, alpha_R, c_R, alpha_rho_R) s_L = 0._wp; s_R = 0._wp diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 37892266a..913838c8d 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -42,16 +42,16 @@ contains end function f_compute_filtered_dtheta !> Computes the mixture coefficients, velocity and pressure of one cell - subroutine s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, vel, vel_sum, qv, j, k, l) + subroutine s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, k, l) $:GPU_ROUTINE(function_name='s_compute_cell_state',parallelism='[seq]', cray_inline=True) type(scalar_field), intent(in), dimension(sys_size) :: q_prim_vf #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), intent(inout), dimension(3) :: alpha + real(wp), intent(inout), dimension(3) :: alpha, alpha_rho real(wp), intent(inout), dimension(3) :: vel #:else - real(wp), intent(inout), dimension(num_fluids) :: alpha + real(wp), intent(inout), dimension(num_fluids) :: alpha, alpha_rho real(wp), intent(inout), dimension(num_vels) :: vel #:endif real(wp), intent(inout) :: rho, gamma, pi_inf, vel_sum, pres @@ -59,9 +59,9 @@ contains integer, intent(in) :: j, k, l real(wp), dimension(2), intent(inout) :: Re #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_rho, Gs + real(wp), dimension(3) :: Gs #:else - real(wp), dimension(num_fluids) :: alpha_rho, Gs + real(wp), dimension(num_fluids) :: Gs #:endif real(wp) :: G_local integer :: i diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 9a04225b0..b680be87e 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -641,11 +641,11 @@ contains real(wp) :: rho !< Cell-avg. density #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: vel !< Cell-avg. velocity - real(wp), dimension(3) :: alpha !< Cell-avg. volume fraction + real(wp), dimension(3) :: vel !< Cell-avg. velocity + real(wp), dimension(3) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density #:else - real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity - real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction + real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity + real(wp), dimension(num_fluids) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density #:endif real(wp) :: vel_sum !< Cell-avg. velocity sum real(wp) :: pres !< Cell-avg. pressure @@ -664,19 +664,21 @@ contains end if dt_local = huge(1.0_wp) - $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, fl, max_dt]', & - & reduction='[[dt_local]]', reductionOp='[min]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, alpha_rho, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, fl, & + & max_dt]', reduction='[[dt_local]]', reductionOp='[min]') do l = 0, p do k = 0, n do j = 0, m if (igr) then - call s_compute_cell_state(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, alpha, vel, vel_sum, qv, j, k, l) + call s_compute_cell_state(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, & + & qv, j, k, l) else - call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, vel, vel_sum, qv, j, k, l) + call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, Re, alpha, alpha_rho, vel, vel_sum, qv, j, & + & k, l) end if ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c, alpha_rho) if (any_non_newtonian) then Re(1) = 0._wp diff --git a/tests/5AC2F65D/golden-metadata.txt b/tests/5AC2F65D/golden-metadata.txt new file mode 100644 index 000000000..3a66558f5 --- /dev/null +++ b/tests/5AC2F65D/golden-metadata.txt @@ -0,0 +1,191 @@ +This file was created on 2026-09-02 20:24:12.370246. + +mfc.sh: + + Invocation: test -j 1 --gpu mp --generate --only 5AC2F65D + Lock: mpi=Yes & gpu=Mp & debug=No & reldebug=No & gcov=No & unified=No & single=No & mixed=No & fastmath=No + Git: 3ffde53b4c166aa8e15348abe433630a7e148ef4 on feature/mie-gruneisen-solver (dirty) + +syscheck: + + CMake Configuration: + + CMake v3.25.2 on k006-004-v7.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : LLVMFlang v23.0.0 (/work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang) + + PRE_PROCESS : OFF + SIMULATION : OFF + POST_PROCESS : OFF + SYSCHECK : ON + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : ON + + Fypp : /home1/sbryngelson/work/mfc-eos/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc + CXX : /opt/ohpc/pub/compiler/gcc/12.2.0/bin/c++ + FC : /work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +simulation: + + CMake Configuration: + + CMake v3.25.2 on k006-007-v8.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : LLVMFlang v23.0.0 (/work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang) + + PRE_PROCESS : OFF + SIMULATION : ON + POST_PROCESS : OFF + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : ON + + Fypp : /home1/sbryngelson/work/mfc-eos/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : + CXX : + FC : /work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +pre_process: + + CMake Configuration: + + CMake v3.25.2 on k006-007-v8.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : LLVMFlang v23.0.0 (/work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang) + + PRE_PROCESS : ON + SIMULATION : OFF + POST_PROCESS : OFF + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : ON + + Fypp : /home1/sbryngelson/work/mfc-eos/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : + CXX : + FC : /work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang + OMPI_CC : + OMPI_CXX : + OMPI_FC : + +post_process: + + CMake Configuration: + + CMake v3.25.2 on k006-007-v8.hpcfund + + C : GNU v12.2.0 (/opt/ohpc/pub/compiler/gcc/12.2.0/bin/cc) + Fortran : LLVMFlang v23.0.0 (/work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang) + + PRE_PROCESS : OFF + SIMULATION : OFF + POST_PROCESS : ON + SYSCHECK : OFF + DOCUMENTATION : OFF + ALL : OFF + + MPI : ON + OpenACC : OFF + OpenMP : ON + + Fypp : /home1/sbryngelson/work/mfc-eos/build/venv/bin/fypp + Doxygen : + + Build Type : Release + + Configuration Environment: + + CC : + CXX : + FC : /work1/spencerbryngelson/sbryngelson/software/therock-afar-23.2.1-gfx90a-7.13.0-7357b5084b/bin/amdflang + 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): 16 + On-line CPU(s) list: 0-15 + Vendor ID: AuthenticAMD + Model name: AMD EPYC 7V13 64-Core Processor + CPU family: 25 + Model: 1 + Thread(s) per core: 1 + Core(s) per socket: 1 + Socket(s): 16 + Stepping: 1 + BogoMIPS: 4890.81 + Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm rep_good nopl cpuid extd_apicid tsc_known_freq pni pclmulqdq ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand hypervisor lahf_lm cmp_legacy svm cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw perfctr_core ssbd ibrs ibpb stibp vmmcall fsgsbase tsc_adjust bmi1 avx2 smep bmi2 invpcid rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves clzero xsaveerptr wbnoinvd arat npt lbrv nrip_save tsc_scale vmcb_clean flushbyasid pausefilter pfthreshold v_vmsave_vmload vgif umip pku ospke vaes vpclmulqdq rdpid overflow_recov succor arch_capabilities + Virtualization: AMD-V + Hypervisor vendor: KVM + Virtualization type: full + L1d cache: 1 MiB (16 instances) + L1i cache: 1 MiB (16 instances) + L2 cache: 8 MiB (16 instances) + L3 cache: 256 MiB (16 instances) + NUMA node(s): 1 + NUMA node0 CPU(s): 0-15 + Vulnerability Gather data sampling: 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 Reg file data sampling: Not affected + Vulnerability Retbleed: Not affected + Vulnerability Spec rstack overflow: Vulnerable: Safe RET, no microcode + Vulnerability Spec store bypass: Mitigation; Speculative Store Bypass disabled via prctl + Vulnerability Spectre v1: Mitigation; usercopy/swapgs barriers and __user pointer sanitization + Vulnerability Spectre v2: Mitigation; Retpolines; IBPB conditional; IBRS_FW; STIBP disabled; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected + Vulnerability Srbds: Not affected + Vulnerability Tsa: Vulnerable: No microcode + Vulnerability Tsx async abort: Not affected + Vulnerability Vmscape: Not affected + diff --git a/tests/5AC2F65D/golden.txt b/tests/5AC2F65D/golden.txt new file mode 100644 index 000000000..9658540c4 --- /dev/null +++ b/tests/5AC2F65D/golden.txt @@ -0,0 +1,24 @@ +D/cons.1.00.000000.dat 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 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 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 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 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.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 +D/cons.1.00.000050.dat 0.81000000000105 0.81000000000473 0.80999999999401 0.80999999996407 0.81000000000449 0.81000000016452 0.80999999981514 0.8099999967751 0.80999998141795 0.80999990100208 0.80999949546846 0.8099975588105 0.80998882699222 0.80995172941667 0.80980352403798 0.80924797225729 0.80730831776682 0.80075942533141 0.78233387677271 0.75588175253221 0.73483290471555 0.71614988442417 0.70241463686298 0.70205934354209 0.69503702803707 0.69342821894165 0.69501935538136 0.69512815632082 0.69361001154088 0.69408521778969 0.6707512794281 0.52516839918889 0.34219691096305 0.29772883960781 0.28876143447198 0.2884605431297 0.28807531149753 0.28815517098807 0.28701819240252 0.28859239274965 0.28413541042879 0.28340274424475 0.27460427844021 0.26210426197247 0.25296423960551 0.25065925679242 0.25014529379589 0.25003097586269 0.25000638261896 0.25000126864968 0.25000024292396 0.25000004482362 0.25000000799052 0.25000000131017 0.25000000008082 0.24999999993363 0.24999999999883 0.25000000001463 0.25000000000223 0.24999999999803 0.24999999999956 0.25000000000028 0.25000000000008 0.24999999999996 0.24999999999999 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 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 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.25000000000001 0.25000000000004 0.24999999999996 0.24999999999973 0.25000000000012 0.25000000000179 0.24999999999994 0.24999999998804 0.24999999998992 0.25000000005015 0.25000000002667 0.24999999932688 0.24999999570048 0.24999997733626 0.24999988618782 0.24999945721195 0.24999755548974 0.24998964114632 0.24995884616444 0.24984741610064 0.24947518055954 0.24834304184915 0.24518456060669 0.23771038458238 0.22630310173052 0.21339612764364 0.20406091387051 0.19532039231634 0.18998694302767 0.1879656297177 0.18508322428614 0.1840754282524 0.18577307269787 0.18484779825405 0.18423423535055 0.18478502416963 0.18343194454611 0.16929362227557 0.13572387472256 0.10967685324411 0.10377754352862 0.10206739439582 0.10185399320593 0.10166403227817 0.10147942942983 0.10131047914388 0.10173444008987 0.10094280213316 0.10206133926594 0.09879577564787 0.09690739495738 0.08962966823484 0.08266732699432 0.08051431750828 0.08010432426208 0.08002065485065 0.0800040168818 0.08000076512497 0.08000014247418 0.08000002590658 0.08000000456683 0.08000000067584 0.08000000002915 0.07999999998024 0.08000000000398 0.08000000000626 0.0800000000004 0.07999999999908 0.07999999999988 0.08000000000015 0.08000000000003 0.07999999999998 0.07999999999999 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 +D/cons.2.00.000000.dat 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 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 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 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 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.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 +D/cons.2.00.000050.dat 0.19000000000025 0.19000000000111 0.1899999999986 0.18999999999157 0.19000000000105 0.19000000003859 0.18999999995664 0.18999999924354 0.18999999564125 0.18999997677827 0.1899998816531 0.1899994273753 0.18999737917101 0.18998867727058 0.18995391304595 0.18982359843072 0.18936861774777 0.18783245779379 0.18351041553929 0.17730559626994 0.1723682122439 0.16798577477913 0.16476392965425 0.16468057091568 0.16303344600589 0.16265568956305 0.16303019460191 0.16305073240852 0.16270144745718 0.16281160689485 0.17013011454079 0.21497356268213 0.26681275625691 0.28701131283212 0.28768718567725 0.28817839050568 0.28830717536174 0.28822609458879 0.28713870024394 0.28857719208347 0.28417355725006 0.28340664506638 0.27460387202306 0.26210418586311 0.25296423976117 0.25065925678732 0.25014529379589 0.25003097586269 0.25000638261896 0.25000126864968 0.25000024292396 0.25000004482362 0.25000000799052 0.25000000131017 0.25000000008082 0.24999999993363 0.24999999999883 0.25000000001463 0.25000000000223 0.24999999999803 0.24999999999956 0.25000000000028 0.25000000000008 0.24999999999996 0.24999999999999 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 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 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.25000000000001 0.25000000000004 0.24999999999996 0.24999999999973 0.25000000000012 0.25000000000179 0.24999999999994 0.24999999998804 0.24999999998992 0.25000000005015 0.25000000002667 0.24999999932688 0.24999999570048 0.24999997733626 0.24999988618782 0.24999945721195 0.24999755548974 0.24998964114632 0.24995884616444 0.24984741610064 0.24947518055954 0.24834304184915 0.24518456060602 0.23771038466488 0.22630310005793 0.21339620634691 0.20406050231277 0.19532650055959 0.18997055155394 0.1880171991449 0.18494202219434 0.18426045859518 0.18556449988766 0.18482398305435 0.18442509204953 0.18461020240768 0.18142631549458 0.15861150695561 0.09815591755868 0.04573431949285 0.03195804958427 0.02912336957079 0.02878325492241 0.02861399174212 0.02854510399969 0.0284742628043 0.02860880963653 0.02837702767502 0.02870478049695 0.02778235117261 0.02725466992345 0.0252086218109 0.02325019063049 0.02264465181162 0.02252934119872 0.02250580917675 0.02250112974801 0.0225002151914 0.02250004007086 0.02250000728623 0.02250000128442 0.02250000019008 0.0225000000082 0.02249999999444 0.02250000000112 0.02250000000176 0.02250000000011 0.02249999999974 0.02249999999997 0.02250000000004 0.02250000000001 0.02249999999999 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 +D/cons.3.00.000000.dat 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 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.000050.dat -1.93e-12 -8.68e-12 1.1e-11 6.589e-11 -7.71e-12 -2.9786e-10 3.5394e-10 5.74758e-09 3.402522e-08 1.8167785e-07 9.2596039e-07 4.48034622e-06 2.05058768e-05 8.858906243e-05 0.00036055040814 0.00137955028644 0.00493166149875 0.01684669184666 0.05000144199825 0.09352349591282 0.13457802036663 0.16017177238021 0.17251106843883 0.19282924791141 0.18964047716135 0.1917434592698 0.19579487615641 0.19436027362994 0.19216276928983 0.19335986867506 0.19099904809838 0.16745060648514 0.13669605232212 0.13165692959926 0.13116610456639 0.13101744337351 0.12947334031707 0.12897396657155 0.12949346531537 0.12375989906637 0.12642820095657 0.10424969204749 0.08349832740146 0.0388974720629 0.00929918003281 0.00205646406425 0.00045263621532 9.647084119e-05 1.987662753e-05 3.95075105e-06 7.5650166e-07 1.3958994e-07 2.491564e-08 4.15483e-09 2.4574e-10 -2.0856e-10 -3.89e-12 4.56e-11 6.93e-12 -6.14e-12 -1.36e-12 8.8e-13 2.6e-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 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 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 -2e-14 -1.2e-13 1.2e-13 8.5e-13 -3.8e-13 -5.59e-12 1.8e-13 3.725e-11 3.124e-11 -1.5628e-10 -9.457e-11 2.15488e-09 1.342214e-08 7.060683e-08 3.5443598e-07 1.69030241e-06 7.61246457e-06 3.225804637e-05 0.00012814689263 0.00047501014823 0.00163249361478 0.00514117019023 0.01484969091201 0.03706485282088 0.06975677777979 0.10235570915869 0.13064120043485 0.14308901390357 0.15613723560419 0.16862906493687 0.16498236245308 0.16825280754807 0.17108839977785 0.16870853399508 0.16763936744753 0.16902838390221 0.1674299785145 0.14994581182637 0.10646471753827 0.07051250354644 0.06205508752626 0.06058140801432 0.06032165666191 0.05947203046611 0.05899307934088 0.05916236971218 0.05847081170925 0.05984762355369 0.0561708829033 0.05465774417312 0.04359544860831 0.02469829578075 0.00644967603782 0.00122126740972 0.00024683519949 4.883367963e-05 9.49555628e-06 1.80863434e-06 3.3679341e-07 6.128884e-08 1.090297e-08 1.81759e-09 6.825e-11 -5.307e-11 8.86e-12 1.489e-11 9.3e-13 -2.17e-12 -2.7e-13 3.5e-13 7e-14 -5e-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 +D/cons.4.00.000000.dat 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 +D/cons.4.00.000050.dat 2.50000000000454 2.50000000002043 2.49999999997413 2.49999999984474 2.50000000001941 2.50000000071091 2.49999999920123 2.49999998606523 2.49999991970717 2.49999957223127 2.49999781992669 2.49998945167825 2.49995172213944 2.49979143370112 2.49915119544517 2.4967528567253 2.48839858170219 2.46041492942671 2.38331724228801 2.27630524710102 2.19554929540433 2.12524215896009 2.0738035741938 2.07685921311746 2.04912756836397 2.04322331502971 2.05138614680461 2.04785573157922 2.04718543777935 2.0482722555776 2.05234243203676 2.08840650846582 2.16323761608239 2.16906680383295 2.17718935317326 2.17552251251586 2.17825807972824 2.17590779026632 2.16616166128837 2.17750363478678 2.1354146786191 2.12281073816574 2.03360971262453 1.90902485812864 1.82188235301559 1.80050641936285 1.79577848372874 1.79472876013775 1.79450302383195 1.79445608787473 1.7944466739483 1.79444485582574 1.79444451777962 1.79444445646891 1.79444444518622 1.79444444383527 1.79444444443368 1.79444444457876 1.79444444446489 1.79444444442636 1.79444444444044 1.79444444444703 1.7944444444452 1.79444444444411 1.79444444444432 1.79444444444449 1.79444444444446 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444444 1.79444444444445 1.79444444444445 1.79444444444443 1.7944444444444 1.79444444444451 1.79444444444479 1.7944444444441 1.79444444444194 1.79444444444557 1.7944444444609 1.79444444444391 1.79444444433468 1.79444444435196 1.79444444490472 1.79444444468923 1.79444443826666 1.79444440498436 1.79444423644166 1.79444339990191 1.79443946286427 1.79442200942756 1.79434937588646 1.79406678495725 1.79304461264412 1.78963406031647 1.7792979034799 1.75075094895182 1.68476375074385 1.58821375391779 1.48464001555593 1.41557655124407 1.35006654561488 1.31295242612592 1.30493834831349 1.27553398330751 1.28036280659909 1.28292288671433 1.27985649858138 1.27717762404129 1.27755677304824 1.26364020321158 1.18574941700837 0.96478004983766 0.78712848759005 0.73644627775421 0.72947187916884 0.7286599688487 0.7276811893064 0.72481956959847 0.72358911686263 0.72587385661699 0.71987513721717 0.72740601099551 0.69817035496179 0.67717114549128 0.60951601131881 0.54926722194352 0.53183442907785 0.52859728666688 0.52793989465871 0.52780930042835 0.52778378193063 0.527778895806 0.52777798107252 0.52777781361472 0.52777778308124 0.52777777800653 0.52777777762269 0.52777777780899 0.52777777782688 0.52777777778092 0.52777777777056 0.52777777777688 0.52777777777893 0.52777777777801 0.52777777777761 0.52777777777773 0.5277777777778 0.52777777777779 0.52777777777777 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 0.52777777777778 +D/cons.5.00.000000.dat 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 +D/cons.5.00.000050.dat 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.87548479251571 0.72967224236376 0.55537713832821 0.50710264386236 0.50078553204384 0.50007097186444 0.50000534443764 0.50000033791608 0.50000001823989 0.50000000088384 0.49999999991387 0.49999999999592 0.50000000000424 0.49999999999994 0.49999999999999 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.49920834538201 0.49127332794032 0.44767194141776 0.32874918238154 0.23234903107108 0.20584022348196 0.20095186797927 0.20013341159197 0.20001627616279 0.20000174730049 0.20000016722647 0.20000001424383 0.20000000109008 0.19999999986005 0.19999999998423 0.2000000000108 0.20000000000013 0.19999999999994 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 +D/cons.6.00.000000.dat 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 +D/cons.6.00.000050.dat 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.12451520748429 0.27032775763624 0.44462286167179 0.49289735613764 0.49921446795616 0.49992902813556 0.49999465556236 0.49999966208392 0.49999998176011 0.49999999911617 0.50000000008613 0.50000000000408 0.49999999999576 0.50000000000007 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.50079165461799 0.50872667205968 0.55232805858224 0.67125081761846 0.76765096892892 0.79415977651804 0.79904813202073 0.79986658840803 0.79998372383721 0.79999825269951 0.79999983277353 0.79999998575617 0.79999999890992 0.80000000013995 0.80000000001577 0.7999999999892 0.79999999999987 0.80000000000006 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 +D/prim.1.00.000000.dat 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 0.81 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 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 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 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.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 +D/prim.1.00.000050.dat 0.81000000000105 0.81000000000473 0.80999999999401 0.80999999996407 0.81000000000449 0.81000000016452 0.80999999981514 0.8099999967751 0.80999998141795 0.80999990100208 0.80999949546846 0.8099975588105 0.80998882699222 0.80995172941667 0.80980352403798 0.80924797225729 0.80730831776682 0.80075942533141 0.78233387677271 0.75588175253221 0.73483290471555 0.71614988442417 0.70241463686298 0.70205934354209 0.69503702803707 0.69342821894165 0.69501935538136 0.69512815632082 0.69361001154088 0.69408521778969 0.6707512794281 0.52516839918889 0.34219691096305 0.29772883960781 0.28876143447198 0.2884605431297 0.28807531149753 0.28815517098807 0.28701819240252 0.28859239274965 0.28413541042879 0.28340274424475 0.27460427844021 0.26210426197247 0.25296423960551 0.25065925679242 0.25014529379589 0.25003097586269 0.25000638261896 0.25000126864968 0.25000024292396 0.25000004482362 0.25000000799052 0.25000000131017 0.25000000008082 0.24999999993363 0.24999999999883 0.25000000001463 0.25000000000223 0.24999999999803 0.24999999999956 0.25000000000028 0.25000000000008 0.24999999999996 0.24999999999999 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 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 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.25000000000001 0.25000000000004 0.24999999999996 0.24999999999973 0.25000000000012 0.25000000000179 0.24999999999994 0.24999999998804 0.24999999998992 0.25000000005015 0.25000000002667 0.24999999932688 0.24999999570048 0.24999997733626 0.24999988618782 0.24999945721195 0.24999755548974 0.24998964114632 0.24995884616444 0.24984741610064 0.24947518055954 0.24834304184915 0.24518456060669 0.23771038458238 0.22630310173052 0.21339612764364 0.20406091387051 0.19532039231634 0.18998694302767 0.1879656297177 0.18508322428614 0.1840754282524 0.18577307269787 0.18484779825405 0.18423423535055 0.18478502416963 0.18343194454611 0.16929362227557 0.13572387472256 0.10967685324411 0.10377754352862 0.10206739439582 0.10185399320593 0.10166403227817 0.10147942942983 0.10131047914388 0.10173444008987 0.10094280213316 0.10206133926594 0.09879577564787 0.09690739495738 0.08962966823484 0.08266732699432 0.08051431750828 0.08010432426208 0.08002065485065 0.0800040168818 0.08000076512497 0.08000014247418 0.08000002590658 0.08000000456683 0.08000000067584 0.08000000002915 0.07999999998024 0.08000000000398 0.08000000000626 0.0800000000004 0.07999999999908 0.07999999999988 0.08000000000015 0.08000000000003 0.07999999999998 0.07999999999999 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 +D/prim.2.00.000000.dat 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 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 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 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 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.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 +D/prim.2.00.000050.dat 0.19000000000025 0.19000000000111 0.1899999999986 0.18999999999157 0.19000000000105 0.19000000003859 0.18999999995664 0.18999999924354 0.18999999564125 0.18999997677827 0.1899998816531 0.1899994273753 0.18999737917101 0.18998867727058 0.18995391304595 0.18982359843072 0.18936861774777 0.18783245779379 0.18351041553929 0.17730559626994 0.1723682122439 0.16798577477913 0.16476392965425 0.16468057091568 0.16303344600589 0.16265568956305 0.16303019460191 0.16305073240852 0.16270144745718 0.16281160689485 0.17013011454079 0.21497356268213 0.26681275625691 0.28701131283212 0.28768718567725 0.28817839050568 0.28830717536174 0.28822609458879 0.28713870024394 0.28857719208347 0.28417355725006 0.28340664506638 0.27460387202306 0.26210418586311 0.25296423976117 0.25065925678732 0.25014529379589 0.25003097586269 0.25000638261896 0.25000126864968 0.25000024292396 0.25000004482362 0.25000000799052 0.25000000131017 0.25000000008082 0.24999999993363 0.24999999999883 0.25000000001463 0.25000000000223 0.24999999999803 0.24999999999956 0.25000000000028 0.25000000000008 0.24999999999996 0.24999999999999 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 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 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.25000000000001 0.25000000000004 0.24999999999996 0.24999999999973 0.25000000000012 0.25000000000179 0.24999999999994 0.24999999998804 0.24999999998992 0.25000000005015 0.25000000002667 0.24999999932688 0.24999999570048 0.24999997733626 0.24999988618782 0.24999945721195 0.24999755548974 0.24998964114632 0.24995884616444 0.24984741610064 0.24947518055954 0.24834304184915 0.24518456060602 0.23771038466488 0.22630310005793 0.21339620634691 0.20406050231277 0.19532650055959 0.18997055155394 0.1880171991449 0.18494202219434 0.18426045859518 0.18556449988766 0.18482398305435 0.18442509204953 0.18461020240768 0.18142631549458 0.15861150695561 0.09815591755868 0.04573431949285 0.03195804958427 0.02912336957079 0.02878325492241 0.02861399174212 0.02854510399969 0.0284742628043 0.02860880963653 0.02837702767502 0.02870478049695 0.02778235117261 0.02725466992345 0.0252086218109 0.02325019063049 0.02264465181162 0.02252934119872 0.02250580917675 0.02250112974801 0.0225002151914 0.02250004007086 0.02250000728623 0.02250000128442 0.02250000019008 0.0225000000082 0.02249999999444 0.02250000000112 0.02250000000176 0.02250000000011 0.02249999999974 0.02249999999997 0.02250000000004 0.02250000000001 0.02249999999999 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 0.0225 +D/prim.3.00.000000.dat 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 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.000050.dat -1.93e-12 -8.68e-12 1.1e-11 6.589e-11 -7.71e-12 -2.9786e-10 3.5394e-10 5.74758e-09 3.402522e-08 1.8167787e-07 9.2596096e-07 4.48035972e-06 2.050615966e-05 8.859434206e-05 0.00036063788552 0.00138083229162 0.00494810436864 0.01704109869222 0.05176967177449 0.10021942114075 0.14834419606721 0.18116198652653 0.19893373187448 0.22247648307745 0.22100804409202 0.22397741315418 0.22818597849066 0.22647990550981 0.22440756487677 0.22565128391769 0.22714148448081 0.22624120116341 0.22445629302752 0.22515459054744 0.22754170967126 0.22720880559958 0.22463094085765 0.22376502199889 0.22553672519457 0.21442553855667 0.22246385003027 0.18392372111935 0.15203402813856 0.07420229914933 0.01838042413514 0.00410211075105 0.00090474661436 0.00019291777919 3.975224016e-05 7.901462e-06 1.51300186e-06 2.7917983e-07 4.983128e-08 8.30966e-09 4.9149e-10 -4.1712e-10 -7.78e-12 9.12e-11 1.387e-11 -1.227e-11 -2.72e-12 1.75e-12 5.1e-13 -2.3e-13 -8e-14 3e-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 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 1e-14 3e-14 -5e-14 -2.3e-13 2.3e-13 1.7e-12 -7.6e-13 -1.117e-11 3.6e-13 7.45e-11 6.249e-11 -3.1256e-10 -1.8914e-10 4.30976e-09 2.684429e-08 1.4121367e-07 7.0887229e-07 3.38061217e-06 1.5225078e-05 6.45187661e-05 0.0002563359821 0.00095060048178 0.0032718557636 0.01035094470928 0.03028267945438 0.07796220783448 0.15412245237505 0.23982555685022 0.32010376141639 0.36628734673953 0.41093342763544 0.44850203783773 0.44586785367302 0.45679178585629 0.4607354935473 0.45637384979171 0.45472704740657 0.45758139721612 0.45889047022211 0.4572841302542 0.45521127113985 0.45371579343132 0.45717623582084 0.46178104450814 0.46174929069734 0.45650086354435 0.45370729496105 0.45584996220748 0.44859102279547 0.46278767643344 0.42955226480033 0.43181034153431 0.35111729697918 0.21507021543869 0.06089338366734 0.01183869340466 0.00240501202394 0.00047630316806 9.263492214e-05 1.764504436e-05 3.28578353e-06 5.9793975e-07 1.0637039e-07 1.773256e-08 6.659e-10 -5.1774e-10 8.643e-11 1.4527e-10 9.05e-12 -2.121e-11 -2.65e-12 3.4e-12 6.9e-13 -4.8e-13 -1.4e-13 7e-14 2e-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 +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 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.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.4.00.000050.dat 1.00000000000287 1.0000000000129 0.99999999998366 0.99999999990196 1.00000000001225 1.00000000044889 0.99999999949563 0.99999999120119 0.99999994930082 0.99999972989458 0.99999862343891 0.99999333947631 0.99996951573307 0.9998683007521 0.99946396667887 0.99794859432097 0.99266108104581 0.97484689636325 0.92495407083625 0.85380612052966 0.7976647647322 0.74826856437013 0.71221448392939 0.71134619832648 0.69304231158453 0.68877300041619 0.6933739520266 0.69220855153039 0.69051478406059 0.69135015017854 0.69158953231122 0.69097983815468 0.69157241682048 0.69200276168445 0.69132316910697 0.69104670771917 0.69190326431733 0.69110430851344 0.6858359234196 0.69270753776087 0.67044900239644 0.66635684320214 0.62183650867798 0.55944592291062 0.51442828918191 0.50319936763442 0.5007046482312 0.50015020481142 0.50003094888565 0.50000615155255 0.50000117791174 0.50000021734478 0.50000003874512 0.50000000635288 0.5000000003919 0.49999999967816 0.49999999999431 0.50000000007096 0.5000000000108 0.49999999999044 0.49999999999789 0.50000000000137 0.5000000000004 0.49999999999982 0.49999999999993 0.50000000000002 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.49999999999999 0.49999999999998 0.50000000000003 0.50000000000018 0.49999999999982 0.49999999999868 0.50000000000059 0.5000000000087 0.49999999999972 0.49999999994201 0.49999999995114 0.50000000024318 0.50000000012933 0.49999999673609 0.49999997915208 0.49999989010603 0.49999944813747 0.49999736808296 0.49998814687755 0.4999497716852 0.49980045880967 0.49926025669749 0.49745653500201 0.49197799170492 0.47675084062861 0.44100822852182 0.38719825181189 0.32746201202618 0.28504510462101 0.24603724814783 0.22228804946476 0.21429842376419 0.19939324099447 0.19941997175649 0.20214765256055 0.20014459168254 0.1984675022006 0.19990803510248 0.20109921961222 0.19988472147928 0.19873164887287 0.19955877467624 0.20025870634696 0.20053956586403 0.2008255522636 0.20047169776999 0.19920181872341 0.19847207680634 0.20003196318441 0.19640901789271 0.20144870577267 0.18601675508783 0.17704289546898 0.14329948046627 0.11176939032515 0.10224808525196 0.10045518016591 0.10009008644388 0.10001751838858 0.10000333680684 0.10000062134616 0.10000011298149 0.10000001991646 0.10000000294741 0.10000000012713 0.09999999991381 0.10000000001735 0.10000000002729 0.10000000000175 0.09999999999599 0.0999999999995 0.10000000000064 0.10000000000013 0.09999999999991 0.09999999999997 0.10000000000001 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 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.5.00.000000.dat 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 +D/prim.5.00.000050.dat 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.87548479251571 0.72967224236376 0.55537713832821 0.50710264386236 0.50078553204384 0.50007097186444 0.50000534443764 0.50000033791608 0.50000001823989 0.50000000088384 0.49999999991387 0.49999999999592 0.50000000000424 0.49999999999994 0.49999999999999 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.49920834538201 0.49127332794032 0.44767194141776 0.32874918238154 0.23234903107108 0.20584022348196 0.20095186797927 0.20013341159197 0.20001627616279 0.20000174730049 0.20000016722647 0.20000001424383 0.20000000109008 0.19999999986005 0.19999999998423 0.2000000000108 0.20000000000013 0.19999999999994 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 +D/prim.6.00.000000.dat 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 +D/prim.6.00.000050.dat 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.12451520748429 0.27032775763624 0.44462286167179 0.49289735613764 0.49921446795616 0.49992902813556 0.49999465556236 0.49999966208392 0.49999998176011 0.49999999911617 0.50000000008613 0.50000000000408 0.49999999999576 0.50000000000007 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.50079165461799 0.50872667205968 0.55232805858224 0.67125081761846 0.76765096892892 0.79415977651804 0.79904813202073 0.79986658840803 0.79998372383721 0.79999825269951 0.79999983277353 0.79999998575617 0.79999999890992 0.80000000013995 0.80000000001577 0.7999999999892 0.79999999999987 0.80000000000006 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 \ No newline at end of file diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 9d7351245..6cc91a9db 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -42,9 +42,14 @@ "check_eos_selector": { "title": "Equation of State Selector", "category": "Thermodynamic Constraints", - "math": r"\Pi_\infty = 0 \;\; \text{for an ideal gas}", + "math": r"\rho e = \Gamma\,p + \Pi(\rho), \quad \Gamma = 1/\Gamma_G, \quad \Pi(\rho) = \rho\, e_{\mathrm{ref}}(\rho) - p_{\mathrm{ref}}(\rho)/\Gamma_G", "explanation": ( - "An ideal gas is the stiffened-gas equation of state with no stiffness, so a case that selects it may not " "set pi_inf at all: the selector determines the stiffness, not the input." + "Every backend supplies the same two coefficients, and a case may only set the parameters its backend reads: an " + "ideal gas has no pi_inf, and a Mie-Gruneisen fluid has neither gamma, pi_inf nor qv. Mie-Gruneisen supplies a linear-Hugoniot " + "reference curve (mg_rho0, mg_c0, mg_s, mg_gruneisen) whose reference energy already carries the formation energy, so " + "qv must be zero; its parameters are read only when that backend is selected." + "The Mie-Gruneisen backend is evaluated per phase along the 5-equation Riemann and time-step paths; " + "the features it is refused with still read stiffened-gas coefficients directly." ), "references": ["Wilfong26"], }, @@ -950,17 +955,63 @@ def check_eos_selector(self): return eos_names = CONSTRAINTS["fluid_pp(1)%eos"]["names"] eos_ideal_gas = eos_names["ideal_gas"] + eos_mg = eos_names["mie_gruneisen"] eos_values = set(eos_names.values()) bub_fac = 1 if self.get("bubbles_euler", "F") == "T" else 0 for i in range(1, num_fluids + 1 + bub_fac): eos = self.get(f"fluid_pp({i})%eos") + mg = {k: self.get(f"fluid_pp({i})%mg_{k}") for k in ("rho0", "c0", "s", "gruneisen")} + # An unset selector is stiffened gas, so stray mg_* parameters must be caught before the early return. + self.prohibit( + (eos if eos is not None else eos_names["stiffened_gas"]) != eos_mg and any(v is not None for v in mg.values()), + f"fluid_pp({i})%mg_* are only read when fluid_pp({i})%eos = 'mie_gruneisen'", + ) if eos is None: continue - self.prohibit(eos not in eos_values, f"fluid_pp({i})%eos must be 'stiffened_gas' or 'ideal_gas'") + self.prohibit(eos not in eos_values, f"fluid_pp({i})%eos must be 'stiffened_gas', 'ideal_gas' or 'mie_gruneisen'") self.prohibit( eos == eos_ideal_gas and self.get(f"fluid_pp({i})%pi_inf") is not None, f"fluid_pp({i})%eos = 'ideal_gas' has no stiffness; do not set fluid_pp({i})%pi_inf", ) + self.prohibit( + eos == eos_mg and any(v is None for v in mg.values()), + f"fluid_pp({i})%eos = 'mie_gruneisen' requires fluid_pp({i})%mg_rho0, mg_c0, mg_s and mg_gruneisen", + ) + if eos == eos_mg and all(v is not None for v in mg.values()): + self.prohibit(mg["rho0"] <= 0 or mg["c0"] <= 0 or mg["gruneisen"] <= 0, f"fluid_pp({i})%mg_rho0, mg_c0 and mg_gruneisen must be positive") + self.prohibit(mg["s"] < 1, f"fluid_pp({i})%mg_s must be >= 1 (u_s = c0 + s u_p; s < 1 gives no shock)") + for k in ("gamma", "pi_inf", "qv"): + self.prohibit( + self.get(f"fluid_pp({i})%{k}") is not None, + f"fluid_pp({i})%{k} is not read with eos = 'mie_gruneisen'; the reference curve replaces it", + ) + # The linear Hugoniot has a pole at mu = 1/(s - 1), its maximum compression. Only the initial + # state can be checked here; the solver does not guard the runtime density. + if mg["s"] > 1: + rho_pole = mg["rho0"] * (1.0 + 1.0 / (mg["s"] - 1.0)) + num_patches = self.get("num_patches", 0) or 0 + for j in range(1, num_patches + 1): + ar = self.get(f"patch_icpp({j})%alpha_rho({i})") + a = self.get(f"patch_icpp({j})%alpha({i})") + if not all(isinstance(v, (int, float)) for v in (ar, a)) or a <= 0: + continue + self.warn( + ar / a > 0.8 * rho_pole, + f"patch_icpp({j}) starts fluid {i} at rho = {ar/a:.4g}, within 20% of the Mie-Gruneisen " + f"Hugoniot pole rho0*s/(s-1) = {rho_pole:.4g}; the reference curve is unphysical beyond it", + ) + if not any(self.get(f"fluid_pp({i})%eos") == eos_mg for i in range(1, num_fluids + 1)): + return + # The per-phase evaluation is wired through the 5-equation Riemann and time-step paths only; every + # feature below still reads the stiffened-gas coefficients directly or mixes without partial densities. + self.prohibit(self.get("model_eqns") != 2, "eos = 'mie_gruneisen' requires model_eqns = 2") + self.prohibit(self.get("riemann_solver") not in (1, 2, 5), "eos = 'mie_gruneisen' requires riemann_solver = 1, 2 or 5") + for flag in ("bubbles_euler", "bubbles_lagrange", "alt_soundspeed", "hypoelasticity", "ib", "igr", "relativity", "mhd", "chemistry", "acoustic_source", "probe_wrt", "c_wrt"): + self.prohibit(self.get(flag, "F") == "T", f"eos = 'mie_gruneisen' is not supported with {flag} = T") + for dir in "xyz": + for bound in ("beg", "end"): + bc = self.get(f"bc_{dir}%{bound}") + self.prohibit(isinstance(bc, int) and -12 <= bc <= -5, f"eos = 'mie_gruneisen' is not supported with characteristic boundary condition bc_{dir}%{bound}") def check_stiffened_eos(self): """Checks constraints on stiffened equation of state fluids parameters""" diff --git a/toolchain/mfc/eos.py b/toolchain/mfc/eos.py new file mode 100644 index 000000000..b77ac254b --- /dev/null +++ b/toolchain/mfc/eos.py @@ -0,0 +1,34 @@ +"""The Mie-Gruneisen reference curve as m_variables_conversion computes it, for tests and validation cases.""" + + +def mg_reference(rho, rho0, c0, s): + """p_ref, e_ref and their d/drho for the linear-Hugoniot curve u_s = c0 + s u_p.""" + mu = rho / rho0 - 1.0 + if mu >= 0.0: + d = 1.0 - (s - 1.0) * mu + p = rho0 * c0**2 * mu * (1.0 + mu) / d**2 + dp_dmu = rho0 * c0**2 * ((1.0 + 2.0 * mu) * d + 2.0 * (s - 1.0) * mu * (1.0 + mu)) / d**3 + else: + p = rho0 * c0**2 * mu + dp_dmu = rho0 * c0**2 + e = p * mu / (2.0 * rho0 * (1.0 + mu)) + de_dmu = (dp_dmu * mu * (1.0 + mu) + p) / (2.0 * rho0 * (1.0 + mu) ** 2) + return p, e, dp_dmu / rho0, de_dmu / rho0 + + +def eos_coefficients(rho, rho0, c0, s, gruneisen): + """Gamma, Pi, dPi/drho: the three numbers s_eos_coefficients returns.""" + p, e, dp, de = mg_reference(rho, rho0, c0, s) + return 1.0 / gruneisen, rho * e - p / gruneisen, e + rho * de - dp / gruneisen + + +def sound_speed(rho, pres, rho0, c0, s, gruneisen): + """c^2 = [((Gamma + 1) p + Pi)/rho - dPi/drho]/Gamma, the frozen single-phase speed the solver uses.""" + gamma, pi, dpi = eos_coefficients(rho, rho0, c0, s, gruneisen) + return ((((gamma + 1.0) * pres + pi) / rho - dpi) / gamma) ** 0.5 + + +def hugoniot_state(u_p, rho0, c0, s): + """Shock speed and density behind a shock of particle velocity u_p driven into the reference state.""" + u_s = c0 + s * u_p + return u_s, rho0 * u_s / (u_s - u_p) diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 15f496211..48a3a0805 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -896,8 +896,8 @@ def _load(): _r(f"{px}sph_har_coeff({ll},{mm})", REAL) # Values must match the hand-written eos_* constants in src/common/m_constants.fpp. - _EOS_NAMES = {"stiffened_gas": 1, "ideal_gas": 2} - _EOS_VALUE_LABELS = {1: "stiffened-gas", 2: "ideal-gas"} + _EOS_NAMES = {"stiffened_gas": 1, "ideal_gas": 2, "mie_gruneisen": 3} + _EOS_VALUE_LABELS = {1: "stiffened-gas", 2: "ideal-gas", 3: "Mie-Gruneisen"} # fluid_pp (10 fluids) # Members present in physical_parameters: gamma, pi_inf, Re, cv, qv, qvp, G. @@ -905,10 +905,12 @@ def _load(): # by upstream #1085/#1093 — they must NOT be registered (namelist read would crash). for f in range(1, NF + 1): px = f"fluid_pp({f})%" - CONSTRAINTS[f"fluid_pp({f})%eos"] = {"choices": [1, 2], "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES} + CONSTRAINTS[f"fluid_pp({f})%eos"] = {"choices": [1, 2, 3], "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES} for a, sym in [("gamma", r"\f$\gamma_k\f$"), ("pi_inf", r"\f$\pi_{\infty,k}\f$"), ("cv", r"\f$c_{v,k}\f$"), ("qv", r"\f$q_{v,k}\f$"), ("qvp", r"\f$q'_{v,k}\f$")]: _r(f"{px}{a}", REAL, math=sym) _r(f"{px}eos", INT, math=r"\f$\mathrm{EOS}_k\f$") + for a, sym in [("mg_rho0", r"\f$\rho_{0,k}\f$"), ("mg_c0", r"\f$c_{0,k}\f$"), ("mg_s", r"\f$s_k\f$"), ("mg_gruneisen", r"\f$\Gamma_{G,k}\f$")]: + _r(f"{px}{a}", REAL, math=sym) _r(f"{px}G", REAL, {"hypoelasticity"}, math=r"\f$G_k\f$") _r(f"{px}Re(1)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (shear)") _r(f"{px}Re(2)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (bulk)") diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index 34a8c4bc5..f0f460fb1 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -7,7 +7,7 @@ from ..state import ARG from .case import CaseGeneratorStack, Nt, TestCaseBuilder, define_case_d, define_case_f, define_convergence_case -from .convergence import ConvergenceSpec, run_amp_sweep, run_dt_sweep, run_h_sweep, run_sod_l1 +from .convergence import ConvergenceSpec, run_amp_sweep, run_dt_sweep, run_h_sweep, run_mg_hugoniot, run_mg_wave_speed, run_sod_l1 # Convergence test specs. # One TestCase per (problem, scheme) pair. Trace prefix "Convergence ->" is @@ -125,6 +125,18 @@ def _h_sweep(case_path, ndim, cons_vars, extra_args, expected, tol, resolutions, ) ) + cases.append( + define_convergence_case( + "Convergence -> Mie-Gruneisen -> acoustic speed", + spec=ConvergenceSpec(runner=run_mg_wave_speed, case_path="examples/1D_mg_acoustic/case.py", expected_order=0.0, tol=1.0e-3, resolutions=[100, 200, 400]), + ) + ) + cases.append( + define_convergence_case( + "Convergence -> Mie-Gruneisen -> Hugoniot", + spec=ConvergenceSpec(runner=run_mg_hugoniot, case_path="examples/1D_mg_impact/case.py", expected_order=0.0, tol=0.005, amps=[0.2, 0.5, 1.0, 1.5]), + ) + ) for label, extra_args, expected, tol, min_N in _CONVERGENCE_SOD_SCHEMES: resolutions = [N for N in _RES_SOD_DEFAULT if min_N is None or N >= min_N] cases.append( @@ -659,6 +671,23 @@ def alter_num_fluids(dimInfo): "patch_icpp(3)%alpha(2)": 0.8, }, ) + if dimInfo[0] == ["x"]: + # Fluid 1 on its own reference curve beside an ideal gas, so one kernel carries both EOS paths. + cases.append( + define_case_d( + stack, + "eos=mie_gruneisen", + { + "fluid_pp(1)%eos": "mie_gruneisen", + "fluid_pp(1)%gamma": None, + "fluid_pp(1)%qv": None, + "fluid_pp(1)%mg_rho0": 0.9, + "fluid_pp(1)%mg_c0": 1.0, + "fluid_pp(1)%mg_s": 1.5, + "fluid_pp(1)%mg_gruneisen": 0.4, + }, + ) + ) if len(dimInfo[0]) > 1: alter_capillary() @@ -2862,6 +2891,8 @@ def foreach_example(): "2D_advection_convergence", "3D_advection_convergence", "2D_hypo_shear_contact", # exercised by the convergence suite + "1D_mg_acoustic", # exercised by the convergence suite + "1D_mg_impact", # exercised by the convergence suite "2D_zero_circ_vortex_analytical", "3D_TaylorGreenVortex_analytical", "3D_IGR_TaylorGreenVortex_nvidia", diff --git a/toolchain/mfc/test/convergence.py b/toolchain/mfc/test/convergence.py index c5a8df8b4..b6fe644b2 100644 --- a/toolchain/mfc/test/convergence.py +++ b/toolchain/mfc/test/convergence.py @@ -38,7 +38,7 @@ import numpy as np -from .. import common +from .. import common, eos CONS_TOL = 1e-10 MFC = ".\\mfc.bat" if os.name == "nt" else "./mfc.sh" @@ -394,6 +394,84 @@ def run_amp_sweep(spec: ConvergenceSpec) -> typing.Tuple[bool, str]: return passed, "\n".join(lines) +def _mg_params(cfg: dict): + return tuple(float(cfg[f"fluid_pp(1)%mg_{k}"]) for k in ("rho0", "c0", "s", "gruneisen")) + + +def run_mg_wave_speed(spec: ConvergenceSpec) -> typing.Tuple[bool, str]: + """Sweep N; the measured speed of a small acoustic pulse must match the analytic Mie-Gruneisen c. + + The derivative term dPi/drho enters the sound speed but not the pressure, so a pulse is the one + observable that isolates it: without that term the speed is off by tens of percent. The residual + is the finite-amplitude correction O(drho/rho) of a simple wave, so the check is absolute, not a rate. + """ + errors = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in spec.resolutions: + cfg, run_dir = _run_mfc(spec.case_path, tmpdir, f"N{N}", ["-N", str(N)] + spec.extra_args, 1) + rho0, c0, s, gruneisen = _mg_params(cfg) + c_exact = eos.sound_speed(rho0, float(cfg["patch_icpp(1)%pres"].split("+")[0]), rho0, c0, s, gruneisen) + Nt = int(cfg["t_step_stop"]) + T = Nt * float(cfg["dt"]) + x_cc = (np.arange(N) + 0.5) / N + peaks = [] + for step in (0, Nt): + rho = _read_field(run_dir, step, 1, 1, N) + i = int(np.argmax(rho)) + if not 0 < i < N - 1: + raise common.MFCException(f"N={N}: pulse peak at the boundary after step {step}") + a, b, c = rho[i - 1], rho[i], rho[i + 1] # parabolic sub-cell peak + peaks.append(x_cc[i] + 0.5 * (a - c) / (a - 2.0 * b + c) / N) + errors.append(abs((peaks[1] - peaks[0]) / T - c_exact) / c_exact) + widths = [6, 16] + lines = [ + f" analytic c = {c_exact:.6f} (need every relative error <= {spec.tol:.1e})", + "", + _table_line(["N", "rel. speed err"], widths), + _table_line(["-" * w for w in widths], widths), + ] + for i, N in enumerate(spec.resolutions): + lines.append(_table_line([str(N), f"{errors[i]:.4e}"], widths)) + lines.append(f"\n Worst relative error: {max(errors):.4e}") + return max(errors) <= spec.tol, "\n".join(lines) + + +def run_mg_hugoniot(spec: ConvergenceSpec) -> typing.Tuple[bool, str]: + """Sweep the closing speed U of a symmetric impact; shock speed and plateau density must sit on the Hugoniot. + + The shock state is set by the jump conditions with the full EOS, so it lands on u_s = c0 + s u_p + only if p_ref and e_ref are the curve the parameters describe. The front position comes from the + integral of the density excess, which is sub-cell accurate. + """ + rows = [] + with tempfile.TemporaryDirectory() as tmpdir: + for U in spec.amps: + tag = f"U{U:.3f}".replace(".", "p") + cfg, run_dir = _run_mfc(spec.case_path, tmpdir, tag, ["--U", str(U)] + spec.extra_args, 1) + rho0, c0, s, gruneisen = _mg_params(cfg) + N, Nt = int(cfg["m"]) + 1, int(cfg["t_step_stop"]) + T = Nt * float(cfg["dt"]) + rho = _read_field(run_dir, Nt, 1, 1, N) + x_cc = (np.arange(N) + 0.5) / N + rho_s = float(np.median(rho[np.abs(x_cc - 0.5) < 0.02])) + x_s = 0.5 + float(np.sum((rho[x_cc > 0.5] - rho0)) / N) / (rho_s - rho0) + u_s, rho_h = eos.hugoniot_state(0.5 * U, rho0, c0, s) + rows.append((U, (x_s - 0.5) / T, u_s - 0.5 * U, rho_s, rho_h)) + widths = [7, 11, 11, 11, 11] + lines = [ + f" (need every relative error <= {spec.tol:.3f})", + "", + _table_line(["U", "D measured", "D Hugoniot", "rho_s meas.", "rho_s Hug."], widths), + _table_line(["-" * w for w in widths], widths), + ] + worst = 0.0 + for U, d_m, d_h, r_m, r_h in rows: + worst = max(worst, abs(d_m - d_h) / d_h, abs(r_m - r_h) / r_h) + lines.append(_table_line([f"{U:.3f}", f"{d_m:.5f}", f"{d_h:.5f}", f"{r_m:.5f}", f"{r_h:.5f}"], widths)) + lines.append(f"\n Worst relative error: {worst:.4e}") + return worst <= spec.tol, "\n".join(lines) + + # Entry point used by test.py. diff --git a/toolchain/mfc/test_case_validator.py b/toolchain/mfc/test_case_validator.py index 4cb7fe6c9..a72d09808 100644 --- a/toolchain/mfc/test_case_validator.py +++ b/toolchain/mfc/test_case_validator.py @@ -389,3 +389,60 @@ def test_not_tripped_without_alt_soundspeed(self): if __name__ == "__main__": unittest.main() + + +class TestMieGruneisenSelector(ConstraintTestCase): + """fluid_pp(i)%eos = 'mie_gruneisen' requires its reference curve and forbids what it makes redundant.""" + + MG = { + "fluid_pp(1)%eos": 3, + "fluid_pp(1)%gamma": None, + "fluid_pp(1)%pi_inf": None, + "fluid_pp(1)%qv": None, + "fluid_pp(1)%mg_rho0": 8930.0, + "fluid_pp(1)%mg_c0": 3940.0, + "fluid_pp(1)%mg_s": 1.49, + "fluid_pp(1)%mg_gruneisen": 2.0, + } + + def warnings_for(self, params): + validator = CaseValidator(dict(params)) + validator.validate("simulation") + return "\n".join(validator.warnings) + + def test_accepts_complete_curve(self): + self.assertAccepts({**BASE, **self.MG}) + + def test_scope_is_the_five_equation_riemann_path(self): + base = {**BASE, **self.MG} + self.assertRejects({**base, "model_eqns": 3}, "requires model_eqns = 2") + self.assertRejects({**base, "alt_soundspeed": "T"}, "not supported with alt_soundspeed = T") + self.assertRejects({**base, "bc_x%beg": -5}, "characteristic boundary condition bc_x%beg") + + def test_requires_all_four_parameters(self): + p = {**BASE, **self.MG} + del p["fluid_pp(1)%mg_s"] + self.assertRejects(p, "requires fluid_pp(1)%mg_rho0, mg_c0, mg_s and mg_gruneisen") + + def test_rejects_parameters_under_stiffened_gas(self): + self.assertRejects({**BASE, "fluid_pp(1)%mg_c0": 3940.0}, "only read when fluid_pp(1)%eos = 'mie_gruneisen'") + + def test_rejects_the_stiffened_gas_parameters(self): + for k in ("gamma", "pi_inf", "qv"): + self.assertRejects({**BASE, **self.MG, f"fluid_pp(1)%{k}": 1.0}, f"fluid_pp(1)%{k} is not read with eos = 'mie_gruneisen'") + + def test_ideal_gas_may_not_set_pi_inf(self): + self.assertRejects({**BASE, "fluid_pp(1)%eos": 2, "fluid_pp(1)%pi_inf": 0.0}, "has no stiffness; do not set fluid_pp(1)%pi_inf") + self.assertAccepts({**BASE, "fluid_pp(1)%eos": 2, "fluid_pp(1)%pi_inf": None, "fluid_pp(1)%qv": None}) + + def test_rejects_slope_below_one(self): + self.assertRejects({**BASE, **self.MG, "fluid_pp(1)%mg_s": 0.9}, "mg_s must be >= 1") + + def test_warns_near_the_hugoniot_pole(self): + # rho_pole = rho0 * s/(s-1) = 8930 * 1.49/0.49 ~ 27157; start at 90% of it + p = {**BASE, **self.MG, "patch_icpp(1)%alpha(1)": 1.0, "patch_icpp(1)%alpha_rho(1)": 0.9 * 8930.0 * 1.49 / 0.49} + self.assertIn("Hugoniot pole", self.warnings_for(p)) + + def test_no_pole_warning_at_reference_density(self): + p = {**BASE, **self.MG, "patch_icpp(1)%alpha(1)": 1.0, "patch_icpp(1)%alpha_rho(1)": 8930.0} + self.assertNotIn("Hugoniot pole", self.warnings_for(p)) diff --git a/toolchain/mfc/test_eos_mie_gruneisen.py b/toolchain/mfc/test_eos_mie_gruneisen.py new file mode 100644 index 000000000..d0374769a --- /dev/null +++ b/toolchain/mfc/test_eos_mie_gruneisen.py @@ -0,0 +1,80 @@ +"""Manufactured checks for the Mie-Gruneisen backend in m_variables_conversion. + +The Fortran is s_eos_coefficients: a linear-Hugoniot reference curve mapped onto MFC's +rho e = Gamma p + Pi form. These tests pin the maths that mapping rests on, without a solver run. +""" + +import pytest + +from mfc.eos import eos_coefficients, mg_reference + +# Copper-like, no calibrated material: order-of-magnitude values only. +RHO0, C0, S, G0 = 8930.0, 3940.0, 1.49, 2.0 +STATES = [(r, p) for r in (7500.0, 8930.0, 9800.0, 11000.0) for p in (1e8, 5e9, 2e10)] + + +def _fd(f, x, h): + return (f(x + h) - f(x - h)) / (2.0 * h) + + +@pytest.mark.parametrize("rho", [7500.0, 8931.0, 9800.0, 11000.0]) +def test_reference_curve_derivatives(rho): + """The analytic dp_ref/drho and de_ref/drho match central differences of the curve itself.""" + _, _, dp, de = mg_reference(rho, RHO0, C0, S) + h = rho * 1e-6 + dp_fd = _fd(lambda r: mg_reference(r, RHO0, C0, S)[0], rho, h) + de_fd = _fd(lambda r: mg_reference(r, RHO0, C0, S)[1], rho, h) + assert dp == pytest.approx(dp_fd, rel=1e-8) + assert de == pytest.approx(de_fd, rel=1e-8) + + +@pytest.mark.parametrize("rho,pres", STATES) +def test_gamma_pi_form_is_the_mie_gruneisen_pressure(rho, pres): + """rho e = Gamma p + Pi inverts to exactly p = p_ref + rho Gamma_G (e - e_ref).""" + gamma, pi, _ = eos_coefficients(rho, RHO0, C0, S, G0) + e = (gamma * pres + pi) / rho # the energy MFC stores for this p + p_ref, e_ref, _, _ = mg_reference(rho, RHO0, C0, S) + p_mg = p_ref + rho * G0 * (e - e_ref) + assert p_mg == pytest.approx(pres, rel=1e-12) + + +@pytest.mark.parametrize("rho,pres", STATES) +def test_sound_speed_matches_isentrope(rho, pres): + """c^2 = [((Gamma+1)p + Pi)/rho - dPi/drho - p dGamma/drho]/Gamma equals (dp/drho)_s, integrated numerically.""" + gamma, pi, dpi = eos_coefficients(rho, RHO0, C0, S, G0) + c2 = (((gamma + 1.0) * pres + pi) / rho - dpi) / gamma # dGamma/drho = 0 for constant Gamma_G + + # Walk the isentrope: de = p/rho^2 drho, then re-evaluate p at the new state. + e = (gamma * pres + pi) / rho + h = rho * 1e-6 + + def p_at(r, e_): + g, pi_, _ = eos_coefficients(r, RHO0, C0, S, G0) + return (r * e_ - pi_) / g + + c2_fd = (p_at(rho + h, e + pres / rho**2 * h) - p_at(rho - h, e - pres / rho**2 * h)) / (2.0 * h) + assert c2 > 0.0 + assert c2 == pytest.approx(c2_fd, rel=1e-6) + + +def test_stiffened_gas_is_the_degenerate_member(): + """p_ref = -gamma pi_inf, e_ref = 0, Gamma_G = gamma - 1 reproduces MFC's stored gammas and pi_infs.""" + gam, pinf = 4.4, 6.0e8 # water-like stiffened gas + Gamma_G = gam - 1.0 + p_ref, e_ref = -gam * pinf, 0.0 + rho = 1000.0 + Gamma = 1.0 / Gamma_G + Pi = rho * e_ref - p_ref / Gamma_G + assert Gamma == pytest.approx(1.0 / (gam - 1.0), rel=1e-15) # what MFC stores as gammas(i) + assert Pi == pytest.approx(gam * pinf / (gam - 1.0), rel=1e-15) # what MFC stores as pi_infs(i) + + +def test_release_branch_is_c1_at_rho0(): + """Compression and release branches meet with matching p_ref and dp_ref at mu = 0.""" + eps = 1e-9 * RHO0 + p_up, _, dp_up, _ = mg_reference(RHO0 + eps, RHO0, C0, S) + p_dn, _, dp_dn, _ = mg_reference(RHO0 - eps, RHO0, C0, S) + assert p_up == pytest.approx(C0**2 * eps, rel=1e-6) # linear through zero from above ... + assert p_dn == pytest.approx(-(C0**2) * eps, rel=1e-6) # ... and from below, same slope + assert dp_up == pytest.approx(dp_dn, rel=1e-6) + assert dp_up == pytest.approx(C0**2, rel=1e-6) # bulk modulus rho0 c0^2 over rho0