Merge branch 'master' into reconnection

sources/equations.F90

@@ -75,9 +75,10 @@ module equations
       real(kind=8), dimension(:,:), intent(out) :: u
       logical     , optional      , intent(in)  :: s
     end subroutine
-    subroutine cons2prim_iface(u, q)
+    subroutine cons2prim_iface(u, q, s)
       real(kind=8), dimension(:,:), intent(in)  :: u
       real(kind=8), dimension(:,:), intent(out) :: q
+      integer                     , intent(out) :: s
     end subroutine
     subroutine fluxspeed_iface(q, u, f, c)
       real(kind=8), dimension(:,:)          , intent(in)  :: q, u
@@ -1412,46 +1413,37 @@ module equations
 !
 !===============================================================================
 !
-  subroutine update_primitive_variables(uu, qq)
+  subroutine update_primitive_variables(uu, qq, status)
 
-! include external procedures and variables
-!
     use coordinates, only : nn => bcells
     use coordinates, only : nb, ne, nbl, neu
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:,:,:), intent(inout) :: uu
     real(kind=8), dimension(:,:,:,:), intent(inout) :: qq
+    integer                         , intent(out)   :: status
 
-! temporary variables
-!
     integer :: i, j, k = 1
-!
+
 !-------------------------------------------------------------------------------
 !
-! update primitive variables
+    status = 0
-!
+
 #if NDIMS == 3
     do k = nb, ne
 #endif /* NDIMS == 3 */
       do j = nb, ne
 
-! convert conserved variables to primitive ones
+        call cons2prim(uu(1:nv,nb:ne,j,k), qq(1:nv,nb:ne,j,k), status)
-!
+
-        call cons2prim(uu(1:nv,nb:ne,j,k), qq(1:nv,nb:ne,j,k))
+        if (status /= 0) go to 100
 
       end do ! j = nb, ne
 #if NDIMS == 3
     end do ! k = nb, ne
 #endif /* NDIMS == 3 */
 
-! fill out the borders
-!
 #if NDIMS == 3
     do i = nbl, 1, -1
       qq(1:nv, i,nb:ne,nb:ne) = qq(1:nv,nb,nb:ne,nb:ne)
@@ -1486,6 +1478,8 @@ module equations
     end do
 #endif /* NDIMS == 3 */
 
+    100 continue
+
 !-------------------------------------------------------------------------------
 !
   end subroutine update_primitive_variables
@@ -1790,54 +1784,51 @@ module equations
 !
 !     u - the input array of conservative variables;
 !     q - the output array of primitive variables;
+!     s - the status flag;
 !
 !===============================================================================
 !
-  subroutine cons2prim_hd_iso(u, q)
+  subroutine cons2prim_hd_iso(u, q, s)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:), intent(in)  :: u
     real(kind=8), dimension(:,:), intent(out) :: q
+    integer                     , intent(out) :: s
 
-! local variables
-!
     integer :: i, p
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
-! start accounting time for variable conversion
-!
     call start_timer(imc)
 #endif /* PROFILE */
 
-! iterate over all positions
+    s = 0
-!
+
     do i = 1, size(u,2)
 
+      if (u(idn,i) > 0.0d+00) then
         q(idn,i) = u(idn,i)
         q(ivx,i) = u(imx,i) / u(idn,i)
         q(ivy,i) = u(imy,i) / u(idn,i)
         q(ivz,i) = u(imz,i) / u(idn,i)
+      else
+        s = 1
+        go to 100
+      end if
 
     end do
 
-! update primitive passive scalars
-!
     if (ns > 0) then
       do p = isl, isu
         q(p,:) = u(p,:) / u(idn,:)
       end do
     end if
 
+    100 continue
+
 #ifdef PROFILE
-! stop accounting time for variable conversion
-!
     call stop_timer(imc)
 #endif /* PROFILE */
 
@@ -2216,59 +2207,59 @@ module equations
 !
 !     u - the input array of conservative variables;
 !     q - the output array of primitive variables;
+!     s - the status flag;
 !
 !===============================================================================
 !
-  subroutine cons2prim_hd_adi(u, q)
+  subroutine cons2prim_hd_adi(u, q, s)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:), intent(in)  :: u
     real(kind=8), dimension(:,:), intent(out) :: q
+    integer                     , intent(out) :: s
 
-! local variables
-!
     integer      :: i, p
     real(kind=8) :: ek, ei
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
-! start accounting time for variable conversion
-!
     call start_timer(imc)
 #endif /* PROFILE */
 
-! iterate over all positions
+    s = 0
-!
+
     do i = 1, size(u,2)
 
+      if (u(idn,i) > 0.0d+00) then
         q(idn,i) = u(idn,i)
         q(ivx,i) = u(imx,i) / u(idn,i)
         q(ivy,i) = u(imy,i) / u(idn,i)
         q(ivz,i) = u(imz,i) / u(idn,i)
-      ek       = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i)          &
+        ek       = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i))
-                                                + u(imz,i) * q(ivz,i))
         ei       = u(ien,i) - ek
+        if (ei > 0.0d+00) then
           q(ipr,i) = gammam1 * ei
-
+        else
+          s = 1
+          go to 100
+        end if
+      else
+        s = 1
+        go to 100
+      end if
     end do
 
-! update primitive passive scalars
-!
     if (ns > 0) then
       do p = isl, isu
         q(p,:) = u(p,:) / u(idn,:)
       end do
     end if
 
+    100 continue
+
 #ifdef PROFILE
-! stop accounting time for variable conversion
-!
     call stop_timer(imc)
 #endif /* PROFILE */
 
@@ -2692,36 +2683,31 @@ module equations
 !
 !     u - the input array of conservative variables;
 !     q - the output array of primitive variables;
+!     s - the status flag;
 !
 !===============================================================================
 !
-  subroutine cons2prim_mhd_iso(u, q)
+  subroutine cons2prim_mhd_iso(u, q, s)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:), intent(in)  :: u
     real(kind=8), dimension(:,:), intent(out) :: q
+    integer                     , intent(out) :: s
 
-! local variables
-!
     integer :: i, p
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
-! start accounting time for variable conversion
-!
     call start_timer(imc)
 #endif /* PROFILE */
 
-! iterate over all positions
+    s = 0
-!
+
     do i = 1, size(u,2)
 
+      if (u(idn,i) > 0.0d+00) then
         q(idn,i) = u(idn,i)
         q(ivx,i) = u(imx,i) / u(idn,i)
         q(ivy,i) = u(imy,i) / u(idn,i)
@@ -2730,20 +2716,22 @@ module equations
         q(iby,i) = u(iby,i)
         q(ibz,i) = u(ibz,i)
         q(ibp,i) = u(ibp,i)
+      else
+        s = 1
+        go to 100
+      end if
 
     end do
 
-! update primitive passive scalars
-!
     if (ns > 0) then
       do p = isl, isu
         q(p,:) = u(p,:) / u(idn,:)
       end do
     end if
 
+    100 continue
+
 #ifdef PROFILE
-! stop accounting time for variable conversion
-!
     call stop_timer(imc)
 #endif /* PROFILE */
 
@@ -3308,37 +3296,32 @@ module equations
 !
 !     u - the input array of conservative variables;
 !     q - the output array of primitive variables;
+!     s - the status flag;
 !
 !===============================================================================
 !
-  subroutine cons2prim_mhd_adi(u, q)
+  subroutine cons2prim_mhd_adi(u, q, s)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:), intent(in)  :: u
     real(kind=8), dimension(:,:), intent(out) :: q
+    integer                     , intent(out) :: s
 
-! local variables
-!
     integer      :: i, p
     real(kind=8) :: ei, ek, em
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
-! start accounting time for variable conversion
-!
     call start_timer(imc)
 #endif /* PROFILE */
 
-! iterate over all positions
+    s = 0
-!
+
     do i = 1, size(u,2)
 
+      if (u(idn,i) > 0.0d+00) then
         q(idn,i) = u(idn,i)
         q(ivx,i) = u(imx,i) / u(idn,i)
         q(ivy,i) = u(imy,i) / u(idn,i)
@@ -3347,25 +3330,31 @@ module equations
         q(iby,i) = u(iby,i)
         q(ibz,i) = u(ibz,i)
         q(ibp,i) = u(ibp,i)
-      ek       = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i)          &
+        ek       = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i))
-                                                + u(imz,i) * q(ivz,i))
         em       = 0.5d+00 * sum(q(ibx:ibz,i) * q(ibx:ibz,i))
         ei       = u(ien,i) - (ek + em)
+        if (ei > 0.0d+00) then
           q(ipr,i) = gammam1 * ei
+        else
+          s = 1
+          go to 100
+        end if
+      else
+        s = 1
+        go to 100
+      end if
 
     end do
 
-! update primitive passive scalars
-!
     if (ns > 0) then
       do p = isl, isu
         q(p,:) = u(p,:) / u(idn,:)
       end do
     end if
 
+    100 continue
+
 #ifdef PROFILE
-! stop accounting time for variable conversion
-!
     call stop_timer(imc)
 #endif /* PROFILE */
 
@@ -3987,45 +3976,35 @@ module equations
 !
 !     u - the input array of conservative variables;
 !     q - the output array of primitive variables;
+!     s - the status flag;
 !
 !===============================================================================
 !
-  subroutine cons2prim_srhd_adi(u, q)
+  subroutine cons2prim_srhd_adi(u, q, s)
 
-! include external procedures
-!
     use iso_fortran_env, only : error_unit
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:), intent(in)  :: u
     real(kind=8), dimension(:,:), intent(out) :: q
+    integer                     , intent(out) :: s
 
-! local variables
-!
     logical      :: info
     integer      :: i, p
     real(kind=8) :: mm, bb, mb, en, dn
     real(kind=8) :: w , vv, vm, vs
 
-! local parameters
-!
     character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srhd_adi()'
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
-! start accounting time for variable conversion
-!
     call start_timer(imc)
 #endif /* PROFILE */
 
-! iterate over all positions
+    s = 0
-!
+
     do i = 1, size(u,2)
 
 ! prepare variables which do not change during the Newton-Ralphson iterations
@@ -4062,9 +4041,8 @@ module equations
         write(error_unit,"(a,5(1x,1es24.16e3))") "U          = ", u(:,i)
         write(error_unit,"(a,3(1x,1es24.16e3))") "D, |m|², E = ", dn, mm, en
 
-! set pressure to zero so we can hopefully fix it later
+        s = 1
-!
+        go to 100
-        q(ipr,i) = 0.0d+00
 
       end if
 
@@ -4078,9 +4056,9 @@ module equations
       end do
     end if
 
+    100 continue
+
 #ifdef PROFILE
-! stop accounting time for variable conversion
-!
     call stop_timer(imc)
 #endif /* PROFILE */
 
@@ -5126,43 +5104,35 @@ module equations
 !
 !     u - the input array of conservative variables;
 !     q - the output array of primitive variables;
+!     s - the status flag;
 !
 !===============================================================================
 !
-  subroutine cons2prim_srmhd_adi(u, q)
+  subroutine cons2prim_srmhd_adi(u, q, s)
 
-! include external procedures
-!
     use iso_fortran_env, only : error_unit
 
-! local variables are not implicit by default
-!
     implicit none
 
-! input/output arguments
-!
     real(kind=8), dimension(:,:), intent(in)  :: u
     real(kind=8), dimension(:,:), intent(out) :: q
+    integer                     , intent(out) :: s
 
-! local variables
-!
     logical      :: info
     integer      :: i, p
     real(kind=8) :: mm, mb, bb, en, dn
     real(kind=8) :: w, wt, vv, vm, vs, fc
 
-! local parameters
-!
     character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srmhd_adi()'
-!
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
-! start accounting time for variable conversion
-!
     call start_timer(imc)
 #endif /* PROFILE */
 
+    s = 0
+
 ! iterate over all positions
 !
     do i = 1, size(u,2)
@@ -5216,9 +5186,8 @@ module equations
           write(error_unit,"(a,6(1x,1es24.16e3))") "D, |m|², m.B, |B|², E, W = "&
                           , dn, mm, mb, bb, en, w
 
-! set pressure to zero so we can hopefully fix it later
+          s = 1
-!
+          go to 100
-          q(ipr,i) = 0.0d+00
 
         end if ! p <= 0
 
@@ -5231,25 +5200,22 @@ module equations
         write(error_unit,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = "     &
                         , dn, mm, mb, bb, en
 
-! set pressure to zero so we can hopefully fix it later
+        s = 1
-!
+        go to 100
-        q(ipr,i) = 0.0d+00
 
       end if ! unphysical state
 
     end do
 
-! update primitive passive scalars
-!
     if (ns > 0) then
       do p = isl, isu
         q(p,:) = u(p,:) / u(idn,:)
       end do
     end if
 
+    100 continue
+
 #ifdef PROFILE
-! stop accounting time for variable conversion
-!
     call stop_timer(imc)
 #endif /* PROFILE */