Return to Revisions

First of all, I'd like to emphasize again that it's important to make sure the equation system itself is correct. Anyway, after the correction the system is less demanding, so let me post an answer.

The idea is simple: the system is nonlinear, derivative of $S_0$ in $\theta$ direction of firsteq is of 1st order, derivative of $S_1$ in $\theta$ direction of secondeq is of 1st order, the solution of the system seems to be stiff, so the experience obtained in e.g. this post may help. Let's try artificial viscosity.

After adding artificial viscosity, something related to the issue discussed in this post seems to show up, so I use pdetoode to discretize the system.

V[ϕ_, θ_] := 
  0.045360921162651446/Sqrt[1.00001 + Cos[ϕ] Sin[θ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 - ϕ]] + 
   0.045360921162651446/Sqrt[1.00001 - Sin[θ] Sin[π/6 + ϕ]];

Clear[points, grid, domain]
points@ϕ = points@θ = 100;
eps = 0(*Pi/10^2*);
domain@ϕ = {0, 2 Pi}; domain@θ = {eps, Pi - eps};
(grid@# = Array[# &, points@#, domain@#]) & /@ {ϕ, θ};

difforder = 2;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = 
  pdetoode[{Subscript[S, 0], Subscript[S, 1]}[ϕ, θ, t], t, 
   grid /@ {ϕ, θ}, difforder, {True, False}];
μ0 = 10^-2; μ1 = 5 10^-2;
Unevaluated[
   eq = {μ0 (D[Subscript[S, 0], {θ, 2}] + 
          D[Subscript[S, 0], {ϕ, 2}]) - 
       D[Subscript[S, 0], 
        t] == -V[ϕ, θ] + (D[Subscript[S, 0], θ]^2 + 
         Csc[θ]^2 D[Subscript[S, 0], ϕ]^2), μ1 (D[Subscript[S, 
           1], {θ, 2}] + D[Subscript[S, 1], {ϕ, 2}]) - 
       D[Subscript[S, 1], 
        t] == -(-Cot[θ] D[Subscript[S, 0], θ] + 
          D[Subscript[S, 0], {θ, 2}] + 
          Csc[θ]^2 D[Subscript[S, 0], {ϕ, 2}]) + 
       2 (D[Subscript[S, 0], θ] D[Subscript[S, 1], θ] + 
          Csc[θ]^2 D[Subscript[S, 0], ϕ] D[Subscript[S, 1], ϕ])};
   ic = {Subscript[S, 0] == Sin[θ]^2, Subscript[S, 1] == Sin[θ]^2} /. 
     t -> 0;
   bc = {D[Subscript[S, 0], θ] == 0, 
      D[Subscript[S, 1], θ] == 
       0} /. {{θ -> domain[θ][[1]]}, {θ -> 
        domain[θ][[-1]]}}] /. {Subscript[S, 0] -> 
    Subscript[S, 0][ϕ, θ, t], 
   Subscript[S, 1] -> Subscript[S, 1][ϕ, θ, t]};
del = #[[2 ;; -2]] &;
ode = Map[del, ptoofunc@eq, {2}] // Quiet;

odebc = ptoofunc@diffbc[t][bc];
odeic = ptoofunc@ic;

varlst = Outer[#@##2 &, {Subscript[S, 0], Subscript[S, 1]}, grid@ϕ, 
   grid@θ];

tend = 10; 

showStatus[status_]:=LinkWrite[$ParentLink,
  SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]
 sollst = NDSolveValue[{ode, odeic, odebc}, varlst, {t, 0, tend}, jianshi[t], 
   Method -> {"EquationSimplification" -> "MassMatrix"}(*,SolveDelayed->
   True*)]; // AbsoluteTiming
(* {251.958, Null} *)

{sols0, sols1} = rebuild[#, grid /@ {ϕ, θ}, -1] & /@ sollst;

The showStatus is from this post.

I've solved $S_0$ and $S_1$ all in one, but since firsteq only involves $S0$ (and further trial suggests it's easier to solve than secondeq), you may consider solve firsteq first.

Let's check the solution. $S_0$:

piclst0 = Table[
   Plot3D[sols0[f, th, t], {th, #3, #4}, {f, #1, #2}, PlotPoints -> 50] & @@ 
    Flatten[domain /@ {ϕ, θ}], {t, 0, tend, tend/25}];

ListAnimate[piclst0]

$S_1$:

piclst1 = Table[
   SphericalPlot3D[sols1[f, th, t], {th, #3, #4}, {f, #1, #2}, PlotPoints -> 50] & @@
     Flatten[domain /@ {ϕ, θ}], {t, 0, tend, tend/25}];

ListAnimate[piclst1]

The solution around $\theta=\pi/2$ looks a bit unusual, I'm not sure if it's the nature of the solution or numeric error. You may adjust points@ϕ, points@θ, μ0, μ1 further to see how the solution changes.