Dear all,

I'm trying to solve a wave equation that describes radial oscillations of a TOV star, depending on radius r and t in Llama multipatch coordinates (Thornburg04).

The equation has a structure of

$\ddot{X} = A X' + B X'' + C X$

I rewrote it as set of coupled equations to evolve it with the MoL thorn:

$\dot{X}=\Pi$

$\dot{H} = \Pi'$

$\dot{\Pi} = AH + B H' +C X$

or as noted in my code as it is attached:

For MoL registered as evolved variables are Xi with rhs Xidot, Pi with rhs Pidot and H with rhs Hdot.

The spatial derivatives drH and drPi are also calculated during the schedule of MoL_CalcRHS in my thorn.

Boundary conditions are applied for r>R and an interval near the TOV radius r=R.

if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND. grid_r(i,j,k) <= (TOV_surface + rprec) )then
          H(i,j,k) = 0
          drPi(i,j,k) = 0
          else if ( grid_r(i,j,k) > (TOV_surface + rprec) ) then
          Xi(i,j,k) = 0
          Pi(i,j,k) = 0
          H(i,j,k) = 0

end if

if (grid_r(i,j,k) > (TOV_surface + rprec) )then
          Xidot(i,j,k) = 0
          Hdot(i,j,k) = 0
          Pidot(i,j,k) = 0
          else if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND. grid_r(i,j,k) <= (TOV_surface + rprec) )then
          Hdot(i,j,k) = 0
          Pidot(i,j,k) = B(i,j,k)*drH(i,j,k)  + C(i,j,k)*Xi(i,j,k)

end if

Problems occur close to the surface of the star. My evolved variables start do diverge close to the surface after a few iterations. Looking at the data it seems that the divergence is  founded by the values of H. I tried following things to encircle the issues:

• If I put Pidot = B*drH + C*Xi the values seem to be ok, whereas for Pidot = A*H the divergence appears. The A-factor only amplifies this behavior. If I put Pidot = H it behaves the same, but much slower.
  
• If I use drXi instead of H (as it is commented out), it does not make any difference.
• If I enlarge the size of the interval (e.g. TOV_surface + 5*rprec) the same divergence appears, but shifted towards grid points next to the interval. Also the divergence appears a few iterations later.
  
• If I change the drXi or H values close to the surface manually (e.g. using a backsided differentiation, or put specific values by hand) it also shifts the divergence (like above).

So far I don't know how to remedy this issue, but maybe I'm overlooking something obvious.

Does anyone have an idea on what I could try?

Thanks a lot!

Best regards and merry Christmas,

Severin Frank