<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Dear all,</p>
<p>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).</p>
<p>The equation has a structure of <br>
</p>
<p><math xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mover><mi>X</mi><mo>̈</mo></mover><mo>=</mo><mi>A</mi><msup><mi>X</mi><mo>′</mo></msup><mo>+</mo><mi>B</mi><msup><mi>X</mi><mo>″</mo></msup><mo>+</mo><mi>C</mi><mi>X</mi></mrow><annotation
encoding="TeX">\ddot{X} = A X' + B X'' + C X</annotation></semantics></math><br>
</p>
<p>I rewrote it as set of coupled equations to evolve it with the
MoL thorn:</p>
<p><math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mover><mi>X</mi><mo>˙</mo></mover><mo>=</mo><mi>Π</mi></mrow></semantics></math></p>
<p><math xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mover><mi>H</mi><mo>˙</mo></mover><mo>=</mo><msup><mi>Π</mi><mo>′</mo></msup></mrow><annotation
encoding="TeX">\dot{H} = \Pi'</annotation></semantics></math></p>
<p><math xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mover><mi>Π</mi><mo>˙</mo></mover><mo>=</mo><mi>A</mi><mi>H</mi><mo>+</mo><mi>B</mi><msup><mi>H</mi><mo>′</mo></msup><mo>+</mo><mi>C</mi><mi>X</mi></mrow><annotation
encoding="TeX">\dot{\Pi} = AH + B H' +C X</annotation></semantics></math><br>
</p>
<p>or as noted in my code as it is attached:<br>
</p>
Pidot = A*H + B*drH + C*Xi<br>
Hdot = drPi<br>
Xidot = Pi <br>
<p>For MoL registered as evolved variables are Xi with rhs Xidot, Pi
with rhs Pidot and H with rhs Hdot. <br>
</p>
<p>The spatial derivatives drH and drPi are also calculated during
the schedule of MoL_CalcRHS in my thorn. <br>
</p>
<p>Boundary conditions are applied for r>R and an interval near
the TOV radius r=R. <br>
</p>
<p><br>
</p>
<p>if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND. grid_r(i,j,k)
<= (TOV_surface + rprec) )then<br>
H(i,j,k) = 0<br>
drPi(i,j,k) = 0<br>
else if ( grid_r(i,j,k) > (TOV_surface + rprec) )
then<br>
Xi(i,j,k) = 0<br>
Pi(i,j,k) = 0<br>
H(i,j,k) = 0</p>
<p>end if<br>
</p>
<p>if (grid_r(i,j,k) > (TOV_surface + rprec) )then<br>
Xidot(i,j,k) = 0<br>
Hdot(i,j,k) = 0<br>
Pidot(i,j,k) = 0<br>
else if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND.
grid_r(i,j,k) <= (TOV_surface + rprec) )then<br>
Hdot(i,j,k) = 0<br>
Pidot(i,j,k) = B(i,j,k)*drH(i,j,k) + C(i,j,k)*Xi(i,j,k)</p>
<p>end if</p>
<p><br>
</p>
<p>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:<br>
</p>
<ul>
<li>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.<br>
</li>
<li>If I use drXi instead of H (as it is commented out), it does
not make any difference.</li>
<li>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.<br>
</li>
<li>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).</li>
</ul>
<p>So far I don't know how to remedy this issue, but maybe I'm
overlooking something obvious.<br>
</p>
<p>Does anyone have an idea on what I could try? <br>
</p>
<p>Thanks a lot!</p>
<p>Best regards and merry Christmas,</p>
<p>Severin Frank<br>
</p>
<math xmlns="http://www.w3.org/1998/Math/MathML"></math>
</body>
</html>