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
Hello Severin,
sorry for the very long delay.
I am no real expert on this others on this list may be better suited.
My worry is that you are assuming that outside of the star (and on the star surface I guess) all time derivatives (and the spatial one in drPi it seems) vanish. This would be true of the TOV solution was actually numerically stationary. For the "usual" Valencia formulations of hydrodynamics this is not quite true though for (at least two reasons): the grid is Cartesian typically leading to interpolation errors when interpolation the correct static TOV data found in the initial data thorn into the grid, this leads to oscillations which can move the stellar surface. Your code may be able to deal with those (and the Tuebingen group are after all experts on this).
The second issue that over time atmosphere actually accretes onto the star thus increasing its mass. This happens b/c while there is a analytic stationary solution the handling of eg the pressure gradient and gravity in the ET (one is a flux, the other is a source term) tends to make it very difficult for the code to find a numerical equilibrium solution so there is going to be some dynamics present in the simulation.
Does you scheme lead to stable solutions when implemented eg in spherical symmetry?
None of this really helps you actually fix the issue I am afraid. The usual solution seems to be to not impose and b/c on the stellar surface but only very far away and accept the fact that the surface will blur out. A sharp surface can (apparently, I have never tried this) be obtained using Wolfgang Kastaun's scheme (which I assume your group is well familiar with) in the Pizza code:
https://publikationen.uni-tuebingen.de/xmlui/bitstream/handle/10900/49026/pd...
Yours, Roland
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
Thank you Roland for this detailed answer!
First of all, some things in my thorn have been changed in the last weeks. Therefore I attached the updated files to this mail. I reformulated the coefficients in my initial data. Now, the code runs much better but is still diverging. However the divergence seems to come from the center of the star instead. Right now I'm using the thornburgnc coordinates and made the spherical part as large as possible as this is the currently most stable configuration for my code.
I attached also a plot of the values of my evolved variable for a point at r=2 for two different time steps, where you can see the behavior of the evolution.
To be honest I'm not 100% sure what you meant with "stable solutions when implemented eg in spherical symmetry?"
What I'm basically trying to do is to have a wave equation that approximates the radial oscillation equations for small amplitudes. The coefficients of this wave equation are filled with static values using the TOVSolver. Then I only evolve the the differential equations, using the MoL thorn. There is no evolution of the GRHydro quantities or the ADM quantities themselves. In principal, time integration of these equations should be possible.
I wonder if I'm using the MoL (or other) thorn(s) wrong for this purpose?
Thank you a lot!
Best regards,
Severin
Hello Severin,
sorry for not responding to your email.
About the stable solution: I was wondering if you had a 1d spherically symmetric toy code where you could have tested that your initial data does indeed lead to almost no evolution or if even in a spherical 1d code you would see material accreting or some oscillation on the surface.
Are you evolution equations for the difference to the static background (TOV solution) valid only in the linear regime or also for non-linear deviations?
My worry would be that even when looking at a difference only, you may be hit by the fact that computing eg a pressure gradient of the TOV solution on the Cactus grid and computing a gravitational "force" will not give perfect balance which would (I expect) show up in your difference equation as driving force potentially shifting where the equilibrium is.
From the look of your plot you seem to be experiencing some sort of instability which causes the amplitude of the oscillations to grow more and more, with smaller timesteps helping a bit but not ultimately curing the issue. You could (his is just a guess) try and play with different MoL integration schemes, eg use the RK2 scheme (very dissipative but stable) to see if this would help (I would stay away from ICN and other non RK schemes as I have no idea how recently those were used for any actual simulations).
In principle you should be able to use MoL for this type of thing. The cases where MoL is difficult to use is trying to use it to evolve a grid scalar or any other quantity that is not defined on the mesh refined grid that Carpet sets up. This affects things like eg integrating this shift at the location of the puncture to follow the puncture (puncture tracking) or trying to integrate particle trajectories (which are basically the same as multiple punctures).
Yours, Roland
Thank you Roland for this detailed answer!
First of all, some things in my thorn have been changed in the last weeks. Therefore I attached the updated files to this mail. I reformulated the coefficients in my initial data. Now, the code runs much better but is still diverging. However the divergence seems to come from the center of the star instead. Right now I'm using the thornburgnc coordinates and made the spherical part as large as possible as this is the currently most stable configuration for my code.
I attached also a plot of the values of my evolved variable for a point at r=2 for two different time steps, where you can see the behavior of the evolution.
To be honest I'm not 100% sure what you meant with "stable solutions when implemented eg in spherical symmetry?"
What I'm basically trying to do is to have a wave equation that approximates the radial oscillation equations for small amplitudes. The coefficients of this wave equation are filled with static values using the TOVSolver. Then I only evolve the the differential equations, using the MoL thorn. There is no evolution of the GRHydro quantities or the ADM quantities themselves. In principal, time integration of these equations should be possible.
I wonder if I'm using the MoL (or other) thorn(s) wrong for this purpose?
Thank you a lot!
Best regards,
Severin
Thank you Roland for your reply!
In the meanwhile I was able to solve the issue. I changed my original evolution requation
ζ̈=Aζ′+Bζ″+Cζ\ddot{\zeta} = A \zeta' + B \zeta'' +C \zeta into a system (1) of
ζ˙=Π\dot{\zeta} = \PiΠ˙=Aζ′+Bζ″+Cζ\dot{\Pi}=A\zeta'+B\zeta''+C\zet I wanted to avoid the second-order numerical derivatives by using a system (2)
ζ˙=Π\dot{\zeta}=\PiΠ˙=AH+BH′+Cζ\dot{\Pi}=AH+BH'+C\zetaH˙=Π′\dot{H}=\Pi'which was source of the problem. This system is analytically equivalent to the one above, but indeed not strongly hyperbolic.
I also had to change the boundaries near the r=0 region of the grid, but the main issue was the aforementioned one. System (1) works with the new boundary conditions.
Thank you very much! Best regards, Severin
users@lists.einsteintoolkit.org
-
Haas, Roland -
Roland Haas
-
Severin Frank