Skip to content

Commit 2cbc654

Browse files
DavidFang03DavidFang03
committed
good units at last??
1 parent 1944714 commit 2cbc654

File tree

2 files changed

+36
-71
lines changed

2 files changed

+36
-71
lines changed

‎src/shammodels/sph/src/modules/ComputeEos.cpp‎

Lines changed: 20 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -419,57 +419,41 @@ void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos_internal
419419
u32 total_elements
420420
= shambase::get_check_ref(storage.part_counts_with_ghost).indexes.get(id);
421421

422-
// UnitSystem<Tscal> unit1{
423-
// si->template get<mega, units::years>(),
424-
// si->template get<units::astronomical_unit>(),
425-
// si->template get<units::kilogramm>() * sol_mass,
426-
// };
427-
// m_e = c = 1
428422
using namespace shamunits;
429423
auto unit_sys = *solver_config.unit_sys;
430-
UnitSystem<Tscal> si{};
431-
Constants<Tscal> ctes_si{si};
432-
Tscal second = si.template get<units::second>(); // une seconde en secondes..
433-
UnitSystem<Tscal> unit2{
434-
si.template get<units::second>() / 1e-3,
435-
si.template get<units::metre>() * (ctes_si.c() * second / 1e-3),
436-
si.template get<units::kilogramm>() * (ctes_si.electron_mass()),
437-
};
438-
Tscal mass
439-
= unit_sys.template to<units::kilogramm>() * unit2.template get<units::kilogramm>();
440-
Tscal length
441-
= unit_sys.template to<units::metre>() * unit2.template get<units::metre>();
442-
Tscal time
443-
= unit_sys.template to<units::second>() * unit2.template get<units::second>();
444-
445-
Tscal betaunit = mass / length / (time * time); // pressure
424+
425+
Tscal mass = unit_sys.template to<units::kilogramm>(); // system mass unit in SI
426+
Tscal length = unit_sys.template to<units::metre>(); // in SI
427+
Tscal time = unit_sys.template to<units::second>(); // in SI
428+
429+
Tscal pressure_unit = mass / length / (time * time);
446430
Tscal density_unit = mass / (length * length * length);
447-
Tscal alphabetaunit = sycl::sqrt(betaunit / sycl::rootn(density_unit, 3));
431+
Tscal velocity_unit = length / time;
432+
433+
// shamlog_debug_ln("CONSTANTS ??", mass, length, time);
434+
// shamlog_debug_ln("RATIOS ??", pressure_unit, density_unit, velocity_unit);
448435

449-
Constants<Tscal> ctes{unit2};
450-
const Tscal m_e = ctes.electron_mass();
451-
const Tscal m_p = ctes.proton_mass();
452-
const Tscal h = ctes.h();
453-
const Tscal c = ctes.c();
454-
shamlog_debug_ln("CONSTANTS ??", "unit2:", "m_e", m_e, "m_p", m_p, "h", h, "c", c);
455-
shamlog_debug_ln("RATIOS ??", betaunit, alphabetaunit);
436+
// static constexpr Tscal pi = shamunits::pi<Tscal>;
437+
// constexpr Tscal coeff_tpf = 2.04399770085e7; // = h / (mp^(1/3) m_e c) SI
438+
// const Tscal alpha = sycl::rootn(3 / (8 * pi), 3) * coeff_tpf;
439+
// constexpr Tscal coeff_p = 5.73180502079e-9; // = m_e^4c^5/h^3 SI
440+
// constexpr Tscal beta = coeff_p * pi / 3;
441+
// shamlog_debug_ln("COEFFS ??", alpha, beta);
456442

457443
sham::kernel_call(
458444
q,
459445
sham::MultiRef{rho_getter},
460446
sham::MultiRef{buf_P, buf_cs},
461447
total_elements,
462-
[mu_e = eos_config->mu_e, unit2, betaunit, alphabetaunit](
448+
[mu_e = eos_config->mu_e, density_unit, pressure_unit, velocity_unit](
463449
u32 i, auto rho, Tscal *__restrict P, Tscal *__restrict cs) {
464450
UnitSystem<Tscal> si{};
465451

466452
using namespace shamrock::sph;
467453
Tscal rho_a = rho(i);
468-
auto const res = EOS::pressure_and_soundspeed(mu_e, rho_a, unit2);
469-
P[i] = res.pressure / betaunit;
470-
cs[i] = res.soundspeed / alphabetaunit;
471-
// coeff_p and coeff_fp are in SI and have respectively
472-
// density**{-1/3} and pressure dimensions
454+
auto const res = EOS::pressure_and_soundspeed(mu_e, rho_a * density_unit);
455+
P[i] = res.pressure / pressure_unit;
456+
cs[i] = res.soundspeed / velocity_unit;
473457
});
474458
});
475459

‎src/shamphys/include/shamphys/eos.hpp‎

Lines changed: 16 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -70,46 +70,27 @@ namespace shamphys {
7070
template<class T>
7171
struct EOS_Fermi {
7272
static constexpr T pi = shamunits::pi<T>;
73-
// static constexpr T h = shamunits::Constants<T>::Si::h;
74-
// static constexpr T m_e = shamunits::Constants<T>::Si::electron_mass;
75-
// static constexpr T m_p = shamunits::Constants<T>::Si::proton_mass;
76-
// static constexpr T c = shamunits::Constants<T>::Si::c;
77-
78-
// hopefully h~1, c=1 = m_e = 1
79-
80-
// static constexpr T coeff_p = pi * 5728e18 / 3.;
81-
// // = pi * m_e * m_e * m_e * m_e * c * c * c * c * c / (3 * h * h * h);
82-
// static constexpr T coeff_pf = 3 * h * h * h / (8 * pi * m_p);
83-
84-
//= \tilde p_F = Fermi momentum divided by m_e*c
85-
// static constexpr T tpf(T mu_e, T rho, T coeff_pf) {
86-
// return sycl::rootn(coeff_pf * rho / mu_e, 3) / (m_e * c);
87-
// }
8873

8974
struct PressureAndCs {
9075
T pressure;
9176
T soundspeed;
9277
};
93-
static constexpr PressureAndCs pressure_and_soundspeed(
94-
T mu_e, T rho, shamunits::UnitSystem<T> unit2) {
95-
96-
shamunits::Constants<T> ctes{unit2};
97-
const T m_e = ctes.electron_mass();
98-
const T m_p = ctes.proton_mass();
99-
const T h = ctes.h();
100-
const T c = ctes.c();
101-
const T coeff_p = pi * sycl::pown(m_e, 4) * sycl::pown(c, 5) / (3 * sycl::pown(h, 3));
102-
const T coeff_pf = h * sycl::rootn(3 / (8 * pi * m_p), 3) / (m_e * c);
103-
104-
// -4*31+8*5+34*31 = 18
105-
// 9.1**4*3**5/6.626**3 = 5728.191760371785
106-
107-
T pf = coeff_pf * sycl::rootn(rho / mu_e, 3);
108-
T pf2 = pf * pf;
109-
110-
T P = coeff_p * (pf * sycl::sqrt(pf2 + 1) * (2 * pf2 - 3) + 3 * sycl::asinh(pf));
111-
T cs2 = 8 * coeff_pf * coeff_p * pf2 * pf2
112-
/ (3 * mu_e * sycl::powr(rho, 2. / 3.) * sycl::sqrt(1 + pf2));
78+
static constexpr PressureAndCs pressure_and_soundspeed(T mu_e, T rho) {
79+
// rho has to be SI !!!!
80+
81+
constexpr T ALPHA
82+
= 0.10064082802851738e-2; //(3/(8pi))**(1./3) * h / (mp^(1/3) m_e c) SI
83+
constexpr T BETA = 6002.332181706928e18; // = (pi/3) * m_e^4c^5/h^3 SI
84+
// constexpr T beta = 6002.332181706928e18;
85+
86+
//\tilde p_F = Fermi momentum divided by m_e*c
87+
const T mu13 = sycl::rootn(mu_e, 3);
88+
T tpf = ALPHA * sycl::rootn(rho, 3) / mu13; // rho has to be SI !!
89+
T tpf2 = tpf * tpf;
90+
91+
T P = BETA * (tpf * sycl::sqrt(tpf2 + 1) * (2 * tpf2 - 3) + 3 * sycl::asinh(tpf));
92+
T cs2 = 8 * ALPHA * BETA * tpf2 * tpf2
93+
/ (3 * mu13 * sycl::powr(rho, 2. / 3.) * sycl::sqrt(1 + tpf2));
11394
return {P, sycl::sqrt(cs2)};
11495
}
11596
};

0 commit comments

Comments
 (0)