1 #include "sph/initial/Stellar.h"
2 #include "objects/geometry/Domain.h"
3 #include "objects/utility/IteratorAdapters.h"
4 #include "physics/Constants.h"
5 #include "physics/Eos.h"
6 #include "quantities/Quantity.h"
7 #include "sph/Materials.h"
8 #include "sph/initial/Distribution.h"
9 #include "system/Factory.h"
10 #include "thread/Scheduler.h"
11
12 NAMESPACE_SPH_BEGIN
13
solveLaneEmden(const Float n,const Float dz,const Float z_max)14 Lut<Float> Stellar::solveLaneEmden(const Float n, const Float dz, const Float z_max) {
15 Float phi = 1._f;
16 Float dphi = 0._f;
17 const Float z0 = 1.e-3_f;
18 Float z = z0;
19 Array<Float> solution;
20 while (phi >= 0._f && z < z_max) {
21 solution.push(phi);
22 const Float d2phi = -2._f / z * dphi - pow(phi, n);
23 dphi += d2phi * dz;
24 phi += dphi * dz;
25 z += dz;
26 }
27 return Lut<Float>(Interval(z0, z), std::move(solution));
28 }
29
polytropicStar(const IEos & eos,const Float radius,const Float mass,const Float n)30 Stellar::Star Stellar::polytropicStar(const IEos& eos, const Float radius, const Float mass, const Float n) {
31 const Float G = Constants::gravity;
32
33 Lut<Float> phi = solveLaneEmden(n);
34 const Float z_star = phi.getRange().upper();
35 const Float dphi_star = phi.derivative()(z_star);
36
37 const Float rho_c = 3._f * mass / (4._f * PI * pow<3>(radius)) * z_star / (-3._f * dphi_star);
38 const Float P_c = G * sqr(mass) / pow<4>(radius) * 1._f / (-4._f * PI * (n + 1) * z_star * dphi_star);
39
40 Array<Float> rho(phi.size()), u(phi.size()), P(phi.size());
41 for (auto el : iterateWithIndex(phi)) {
42 const Size i = el.index();
43 const Float x = el.value().y;
44 rho[i] = rho_c * pow(x, n);
45 P[i] = P_c * pow(x, n + 1);
46 u[i] = eos.getInternalEnergy(rho[i], P[i]);
47 }
48
49 Star star;
50 const Interval range(0._f, radius);
51 star.p = Lut<Float>(range, std::move(P));
52 star.rho = Lut<Float>(range, std::move(rho));
53 star.u = Lut<Float>(range, std::move(u));
54
55 return star;
56 }
57
generateIc(const SharedPtr<IScheduler> & scheduler,const SharedPtr<IMaterial> & material,const IDistribution & distribution,const Size N,const Float radius,const Float mass)58 Storage Stellar::generateIc(const SharedPtr<IScheduler>& scheduler,
59 const SharedPtr<IMaterial>& material,
60 const IDistribution& distribution,
61 const Size N,
62 const Float radius,
63 const Float mass) {
64 SphericalDomain domain(Vector(0._f), radius);
65 Array<Vector> points = distribution.generate(*scheduler, N, domain);
66
67 const RawPtr<EosMaterial> eos = dynamicCast<EosMaterial>(material.get());
68 if (!eos) {
69 throw Exception("Cannot generate IC without equation of state");
70 }
71 const Float gamma = material->getParam<Float>(BodySettingsId::ADIABATIC_INDEX);
72 const Float n = 1._f / (gamma - 1._f);
73 Star star = polytropicStar(eos->getEos(), radius, mass, n);
74
75 const Float rho_min = material->getParam<Interval>(BodySettingsId::DENSITY_RANGE).lower();
76 Array<Float> m(points.size());
77 Array<Float> rho(points.size());
78 Array<Float> u(points.size());
79 Array<Float> p(points.size());
80 const Float v = domain.getVolume() / points.size();
81 for (Size i = 0; i < points.size(); ++i) {
82 const Float r = getLength(points[i]);
83 rho[i] = max(star.rho(r), rho_min);
84 u[i] = star.u(r);
85 p[i] = star.p(r);
86 m[i] = rho[i] * v;
87 }
88
89 /// \todo add to params
90 const Float eta = material->getParam<Float>(BodySettingsId::SMOOTHING_LENGTH_ETA);
91 for (Size i = 0; i < points.size(); ++i) {
92 points[i][H] *= eta;
93 }
94
95 Storage storage(material);
96
97 storage.insert<Vector>(QuantityId::POSITION, OrderEnum::SECOND, std::move(points));
98 storage.insert<Float>(QuantityId::MASS, OrderEnum::ZERO, std::move(m));
99 storage.insert<Float>(QuantityId::ENERGY, OrderEnum::ZERO, std::move(u));
100 storage.insert<Float>(QuantityId::DENSITY, OrderEnum::ZERO, std::move(rho));
101 storage.insert<Float>(QuantityId::PRESSURE, OrderEnum::ZERO, std::move(p));
102 storage.insert<Float>(QuantityId::SOUND_SPEED, OrderEnum::ZERO, 100._f);
103 storage.insert<Size>(QuantityId::FLAG, OrderEnum::ZERO, 0);
104
105 MaterialInitialContext context;
106 context.scheduler = scheduler;
107 context.rng = makeRng<UniformRng>();
108 material->create(storage, context);
109
110 return storage;
111 }
112
113 NAMESPACE_SPH_END
114