Since the borders are undefined, it is better to fill them with the copy
of the last available boundary layers. Even though they will be later
updated by the boundary update, it eliminates the possibilty of
propagating NaNs inside the blocks.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This should not be normally required, but it helps to avoid the
unexpected appearance of NaNs and the code execution termination due
to the floating point exceptions.
This should also help with more advanced open boundary conditions (such
as the conditions which keep the divergence or curl of vector fields
zero).
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
We do not update all elements of array du(:,:,:,:) in
update_increment(), therefore some of its elements may be NaNs. Make
sure that all elements have arithmetic representation by initializing
du(:,:,:,:) with zeros.
Also, indicate that array du(:,:,:,:) is only output argument.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
In nr_initial_brackets_srmhd_adi() we estimate the lower bracket, which
guarantees the physical solution, from the pressure positivity derived
from both, the equation of state and the total energy. We take the
maximum of them as the lower bracket.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
In cons2prim_srmhd_adi() we use the energy equation to calculate the
pressure. If the enthalpy was calculated exactly, this the energy
equation gives pressure equal to the one obtained from the equation of
state. If the physical state couldn't be solved, the returned enthalpy
guarantees that the state is physical (v<1, p>0), but the solution
violates the equation of state.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
The equation for total energy gives another condition for the positivity
of pressure which is a cubic function. This condition is somewhat
stronger than the condition coming from the pressure equation and
costs less computationally.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
If the pressure is negative or zero, print a warning and set it to the
minimum allowed value pmin.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This rewrite removes the dependency on |v|² and its derivative in the
energy function and replaces it with the Lorentz factor. It also
slightly simplifies the energy function.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
If the energy function is positive for the lower bracket, it means that
the solution for energy function is unphysical (p<0). In this case use
the lower bracket which guarantees that the cell is physical.
The total energy should be corrected in this case, but this will be done
later.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
The new subroutine finds the initial brackets and estimates the root
value of the enthalpy which are used in the primitive variables solver
for SRMHD.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
The fluxes of the momenta had the second terms wrongly calculated. The
new calculation uses equations from Mignone & Bodo (2006).
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This change implements the initial bracketing the enthalpy. The minimum
enthalpy is found by solving the pressure positivity equation, which
guarantees also the physicality of velocity. The upper bracket is found
from the energy equation. Then, the brackets are shifted up until the
root region is found.
The 1D and 2D Newton-Raphson methods have been slightly improved as
well. The 1Dw method use the information about brackets to limit the
next step. If the solution shoots out of the brackets, the bisection
method is applied. In the 2Dwv method, we print warnings when the
variables become unphysical.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
If the primitive variable solvers cannot find the physical solution,
print an error message together with the values of the conservative
variables for a problematic cell.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
