NDSolve unable to solve system of PDEs

Question

I'm trying to solve a set of PDEs and am running into problems. The system is defined by:

As can be seen, there are two dynamical equations for $\rho$ and $\phi$, but the equation for $v$ has no time derivative on it. I take Dirichlet boundary conditions for all variables, $\rho(x=-L) = \rho(x=L) = \rho_0$, $\phi(x=-L) = 1, \phi(x=L) = 0$ and $V(x=-L)=v(x=L)=0$. Initial conditions are $\rho(x, t=0)= \rho_0$, $\phi(x, t=0)$ is a sigmoid function going from 1 to 0 and $V(x, t=0)=0$.

Somehow, NDSolve is unable to solve this system under certain conditions. In that case, when I run NDSolve it seems to generate an infinite loop and never terminates, while using up a lot of CPU/memory on my PC.

For reference, here is the current working code thanks to @Alex Trounev (sorry for the poor formatting):

kA := 1/τ ((ρdA - ρ[x, t])/ρdA);
kB := 1/τ ((ρdB - ρ[x, t])/ρdB); 
e := (eM (1 - ϕ[x, t]) + eO ϕ[x, t]) ; 
kD := α e; params = {ρdA -> 0.3, ρdB -> 
   1, ρh -> 1, τ -> 0.5, 
  D1 -> 0.01, η -> 0.1, ξ -> 0.1, kF -> 1, z -> 1, 
  e0 -> 1, α -> 0, eM -> 1, eO -> 1, 
  fc -> 0, χ -> 1, κ -> 
   1}; tmax = 1; L = 1; iv = {ρ[x, 
     0] == (Tanh[8 x] + 1)/2 (ρdA - ρdB) + ρdB, ϕ[
     x, 0] == (1 - Tanh[8 x])/2, v[x, 0] == 0} /. 
  params; bc = {ρ[-L, 
     t] == (Tanh[-8 L] + 1)/
       2 (ρdA - ρdB) + ρdB, ρ[L, 
     t] == (Tanh[8 L] + 1)/
       2 (ρdA - ρdB) + ρdB, ϕ[-L, 
     t] == (1 + Tanh[8 L])/2, ϕ[L, t] == (1 - Tanh[8 L])/2, 
   v[-L, t] == 0, v[L, t] == 0} /. params;
PDEsys = {D[ρ[x, t], t] + 
     D[ρ[x, t] (v[x, t]), 
      x] == (kA (1 - ϕ[x, t]) + kB ϕ[x, t]) ρ[x, t], 
   D[ϕ[x, t], t] + (v[x, t]) D[ϕ[x, t], x] == 
    D1 D [ϕ[x, t], {x, 2}] + 
     z D1 D[ρ[x, t], 
       x] D[ϕ[x, t], x]/ρ[x, t] + (kB - kA + kF) ϕ[x, 
       t] (1 - ϕ[x, t]) + 
     kD (1 - ϕ[x, t]), η D[v[x, t], {x, 2}] - ξ v[x, 
       t] == e D[ρ[x, t], x]/ρ[x, t] + 
     D[e, x] Log[ρ[x, t]/ρh]};
vars = {ρ, ϕ, v}
fullsys = Simplify@Join[PDEsys /. params, bc, iv];
ndsol = NDSolve[fullsys, vars, {x, -L, L}, {t, 0, tmax}, 
  Method -> {"EquationSimplification" -> "Residual", 
    "IndexReduction" -> Automatic}]

Now I'm simplifying a few conditions, and the question is how to get it working when eM is not equal to eO, and when alpha is not equal to zero.

Answer 1

There are exist range of parameters where NDSolve getting solution, for example,

params = {kA -> .1, kB -> 1, D1 -> .1, \[Eta] -> 1, \[Xi] -> 1, 
  kD -> 1, kF -> .5, 
  z -> 1}; rh0 = 1; tmax = 1; L = 1; iv = {\[Rho][x, 0] == 
   rh0, \[Phi][x, 0] == (1 - Tanh[8 x])/2, 
  v[x, 0] == 0}; bc = {\[Rho][-L, t] == rh0, \[Rho][L, t] == 
   rh0, \[Phi][-L, t] == (1 + Tanh[8 L])/2, \[Phi][L, 
    t] == (1 - Tanh[8 L])/2, v[-L, t] == 0, v[L, t] == 0};

PDEsys = {D[\[Rho][x, t], t] + 
     D[\[Rho][x, t] (v[x, t]), 
      x] == (kA (1 - \[Phi][x, t]) + kB \[Phi][x, t]) \[Rho][x, t], 
   D[\[Phi][x, t], t] + (v[x, t]) D[\[Phi][x, t], x] == 
    D1 D[\[Phi][x, t], {x, 2}] + 
     z D1 D[\[Rho][x, t], 
       x] D[\[Phi][x, t], x]/\[Rho][x, t] + (kB - kA + kF) \[Phi][x, 
       t] (1 - \[Phi][x, t]) + 
     kD, \[Eta] D[v[x, t], {x, 2}] - \[Xi] v[x, t] == 
    D[\[Rho][x, t], x]/\[Rho][x, t]};
vars = {\[Rho], \[Phi], v}
fullsys = Join[PDEsys /. params, bc, iv];

ndsol = NDSolve[fullsys, vars, {x, -L, L}, {t, 0, tmax}, 
  Method -> {"EquationSimplification" -> "Residual", 
    "IndexReduction" -> Automatic}]

Visualization

{Plot3D[\[Rho][x, t] /. ndsol[[1]], {x, -L, L}, {t, 0, tmax}, 
  Mesh -> None, ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "\[Rho]"], 
 Plot3D[\[Phi][x, t] /. ndsol[[1]], {x, -L, L}, {t, 0, tmax}, 
  Mesh -> None, ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "\[Phi]"], 
 Plot3D[v[x, t] /. ndsol[[1]], {x, -L, L}, {t, 0, tmax}, Mesh -> None,
   ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "v"]}