I would like to solve relativistic hydrodynamic equations (nonlinear PDEs) introduced here:
I use eqs (33 - 35), (38 - 41), where (40) P(rho)=k*rho^g0 (all with one spatial coordinate "mu" and one temporal "t"). CODE EDITED: 16.04.2020 - this I use.
(*Initial functions-stationary,homogeneous perfect fluid sphere \
structure*)
(****************************************************************)
ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
c = 2.99792*10^10;(*m/s*)
gr = 6.674*10^-8;(*grav. const. in cm^3*g^-1*s^-2*)
gcc = gr/c^2;
m0 = 1.672621*10^-24*gr/c^2;(*proton mass in g trnasformed to cm*)
Ms0 = 1.98855*10^33;
Ms = Ms0*gr/c^2;(*mass of central object in g trnasfomred to cm*)
dr = 10^-5;(*small step and initial m is only e*)
(*initital data*)
g0 = 5/3; rho0 = 10^11; ep0 = 3.64*10^18; e0 =
rho0 (1 + ep0/c^2); pc = (g0 - 1)*rho0*ep0;
dmu = 4*\[Pi]*rho0*dr^2; mumax = 21 Ms0; \[Gamma] = g0; k = pc/rho0^g0;
{pc // N, rho0 // N, e0, ep0 // N, ep0/c^2}
(*Solution TOV and mass equation*)
s = NDSolve[{r'[mu] == Sqrt[1 - 2 m[mu]*gr/(r[mu]*c^2)]/(
4 \[Pi]*rho0*r[mu]^2),
m'[mu] == e0/rho0 Sqrt[1 - (2 m[mu] gcc)/r[mu]], r[dmu] == dr,
m[dmu] == dmu}, {r, m}, {mu, dmu, mumax}];
(*Initial functions to hydrodynamical calculations*)
r0 = r /. s[[1, 1]]; fm0 = m /. s[[1, 2]];
{r0[mumax], fm0[mumax]/Ms0, dmu // N, mumax // N}
f3 = Plot[{fm0[mu]}/Ms0, {mu, dmu, mumax}, Frame -> True,
FrameLabel -> {"\[Mu] [g]", "M/Ms []"}, PlotRange -> All]
f4 = Show[
Plot[{r0[mu]}, {mu, dmu, mumax}, Frame -> True,
FrameLabel -> {"\[Mu] [g]", "r [cm]"}]]
frho0[x_] = If[x < mumax, rho0, 1];
(*Relativistic hydrodynamical equations-collapse of star*)
(**************************************************)
(*introducing of equation*)
G[mu_, t_] = 4 \[Pi]*rho[mu, t]*r[mu, t]^2*D[r[mu, t], mu];(*MW39*)
w[mu_, t_] = 1 + ep[mu, t]/c^2 + p[mu, t]/(rho[mu, t]*c^2);(*MW41*)
a[mu_, t_] = 1/w[mu, t];
ep[mu_, t_] = k*rho[mu, t]^(\[Gamma] - 1)/(\[Gamma] - 1);
p[mu_, t_] = (\[Gamma] - 1) ep[mu, t]*rho[mu, t];(*MW40*)
equt[mu_,
t_] = -a[mu,
t] (4 \[Pi]*r[mu, t]^2*G[mu, t]/w[mu, t]*D[p[mu, t], mu] + (
m[mu, t]*gr)/
r[mu, t]^2 + (4 \[Pi]*gr)/c^2 p[mu, t]*r[mu, t]);(*MW33*)
eqrt[mu_, t_] = a[mu, t]*u[mu, t];(*MW34*)
eqmm[mu_, t_] =
4 \[Pi]*rho[mu, t]*(1 + ep[mu, t]/c^2)*
r[mu, t]^2 D[r[mu, t], mu];(*MW38*)
eqrhort[mu_, t_] = -a[mu, t]*rho[mu, t]*r[mu, t]^2 D[u[mu, t], mu]/
D[r[mu, t], mu];(*MW35*)
(*preparation for solution*)
(*equations*)
eqs = {D[u[mu, t], t] == equt[mu, t], D[r[mu, t], t] == eqrt[mu, t],
D[m[mu, t], mu] == eqmm[mu, t],
D[rho[mu, t]*r[mu, t]^2, t] == eqrhort[mu, t]};
(*boundary conditions*)
bcon = {DirichletCondition[u[mu, t] == 0., mu == dmu],
DirichletCondition[r[mu, t] == r0[dmu], mu == dmu],
DirichletCondition[m[mu, t] == fm0[dmu], mu == dmu],
DirichletCondition[rho[mu, t] == frho0[mumax], mu == mumax]};
(*initial conditions*)
incon = {u[mu, 0] == 0., r[mu, 0] == r0[mu], m[mu, 0] == fm0[mu],
rho[mu, 0] == frho0[mu]};
(*PDEs solution*)
Clear[fu, fr, fm, fro]
{fu, fr, fm, fro} =
NDSolveValue[{eqs, incon, bcon}, {u, r, m, rho}, {mu, dmu,
mumax}, {t, 0, 0.1}]
Initial functions r0[mu], fm0[mu] and frho0[mu] are interpolated functions coming from numerical solution of stationary problem. Result of this solution are error messages:
NDSolveValue::femcnsd: The PDE coefficient -((6.674*10^-8 m[mu])/r[mu]^2)-1.15712*10^-17 r[mu] rho[mu]^(5/3)-3.26355*10^23 r[mu]^4 rho[mu]^(2/3) (r^\[Prime])[mu] (rho^\[Prime])[mu] does not evaluate to a numeric scalar at the coordinate {2.08798*10^34}; it evaluated to Indeterminate instead.
NDSolveValue::femcnsd: The PDE coefficient -((6.674*10^-8 m[mu])/r[mu]^2)-1.15712*10^-17 r[mu] rho[mu]^(5/3)-3.26355*10^23 r[mu]^4 rho[mu]^(2/3) (r^\[Prime])[mu] (rho^\[Prime])[mu] does not evaluate to a numeric scalar at the coordinate {2.08798*10^34}; it evaluated to Indeterminate instead.
Unfortunately, I don't know where is the problem (whole concept, method or...). The problem appears ever in half value of endpoint of integration (mumax/2), doesn't matter what "mumax" is. I'm able to draw (and evaluate in all point of the range) all defined functions in the initial time without problems.
Thank you for help or suggestions.
PS: I'm new here if something is misspelt, marked or unlisted. Please notify me. Thank you.
dmu = 4*\[Pi]*rho0*dr^2. $\endgroup$