If the parallel component of magnetic field is too small, the
intermediate states might produce numerical instabilities, since they
are obtained by the division by a small factor.
In order to remove this situation, we apply full HLLD solver for strong
enough Alfvén wave. If it is weak, the HLLC solver is applied.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
If the parallel component of magnetic field is too small, the
intermediate states might produce numerical instabilities, since they
are obtained by the division by a small factor.
In order to remove this situation, we apply full HLLD solver for strong
enough Alfvén wave. If it is weak, the HLL solver is applied.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
If the parallel component of magnetic field is too small, the
intermediate states might produce numerical instabilities, since they
are obtained by the division by a small factor.
In order to remove this situation, we apply full HLLD solver for strong
enough Alfvén wave. If it is weak, the HLL solver is applied.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
For one of the degenerate situations the state vector was not
calculated. This introduced numerical instabilities. This change fixes
it.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
We should take half of the TVD limited derivatives in order to compare
properly with the high order interpolated derivative. Fix it.
Also, TVD limit both states if the high order interpolation of any of
them overshoot.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
In the sum of the J² we used Fortran subroutine sum() and the sum was
done along the 2nd index. However, the interpretation of this index is
uncertain when the first array rank is 1.
Therefore, instead of using subroutine sum(), the calculation is done
directly.
After this change there are no strange effects appearing on the
boundaries of the blocks.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This additional multi-dimensional limiting can be applied to any
reconstruction method. It is controlled by the parameter 'mlp_limiting'
by setting 'on' or 'off'.
The implementation is based on Gerlinger (2012).
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
In this way the user can define his problem, i.e. the initial setup, the
shape update, the gravitational acceleration profile, in one module.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
If specific boundary condition depends on the position, it would be good
to pass the coordinates to them in order to update the data block.
Also, if boundary conditions are time dependent, we should pass time and
the time increment too.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
If we have time dependent source terms, we should pass tm and dtm to
subroutine update_sources().
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
