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
Ẍ=AX′+BX″+CX\ddot{X} = A X' + B X'' + C X
I rewrote it as set of coupled equations to evolve it with the MoL thorn:
X˙=Π
H˙=Π′\dot{H} = \Pi'
Π˙=AH+BH′+CX\dot{\Pi} = AH + B H' +C X
or as noted in my code as it is attached:
Pidot = A*H + B*drH + C*Xi Hdot = drPi Xidot = Pi
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