Also, calculate vorticity and current density assuming dx=dy=dz=1, so
the criterion value is level independent.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
The previous commit fixes the problem with updating boundaries for
changed blocks. It means that we don't need to mark all their neighbors
for update as well.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
We can have a situation when the whole block was changed and its update
flag was set. The boundaries of the block have to be updated too,
however, the neighbor's update flag is unset. This creates a undesired
situation of the block boundaries not being updated.
Just check if either the block or its neighbor has been set to update to
avoid the above situation.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
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>
