* tweak the way the atmosphere is handled in the general EOS case.

  This is not perfect yet and will need more work.



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@59 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 74 insertions, 34 deletions

diff --git a/src/Whisky_Con2Prim.F90 b/src/Whisky_Con2Prim.F90
index 6179048..5d233f7 100644
--- a/src/Whisky_Con2Prim.F90
+++ b/src/Whisky_Con2Prim.F90
@@ -114,11 +114,14 @@ end module Con2Prim_fortran_interfaces
    @history 
    Trimmed and altered from the GR3D routines, original author Mark Miller.
    2007?: Bruno escluded the points in the atmosphere and excision region from the computation.
-   Aug. 2008: Luca added a check on whether a failure at a given point may be disregarded, because that point will then be restricted from a finer level.
-   This should be completely safe only if *regridding happens at times when all levels are evolved.*
-   Feb. 2009: The above procedure proved to be wrong, so Luca implemented another one. When a failure occurs, it is temporarily ignored, except for storing 
-   the location of where it occured in a mask. At the end, after a Carpet restriction, the mask is checked and if it still contains failures, the run is 
-   aborted with an error message. Only used routines have been updated to use this procedure.
+   Aug. 2008: Luca added a check on whether a failure at a given point may be disregarded, 
+   because that point will then be restricted from a finer level. This should be completely 
+   safe only if *regridding happens at times when all levels are evolved.*
+   Feb. 2009: The above procedure proved to be wrong, so Luca implemented another one. 
+   When a failure occurs, it is temporarily ignored, except for storing the location of where 
+   it occured in a mask. At the end, after a Carpet restriction, the mask is checked and if 
+   it still contains failures, the run is aborted with an error message. Only used routines 
+   have been updated to use this procedure.
    @endhistory 
 
 @@*/
@@ -518,7 +521,8 @@ end subroutine Con2Prim_pt
         It computes the primitive variables on cell boundaries.
         Since reconstruction on conservative had not proved to be very successful, 
         some of the improvements to the C2P routines (e.g. the check about 
-        whether a failure happens in a point that will be restriced anyway) are not implemented here yet.
+        whether a failure happens in a point that will be restriced anyway) 
+        are not implemented here yet.
    @enddesc 
    @calls     
    @calledby   
@@ -1375,6 +1379,11 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
   CCTK_REAL :: psi4pt
   character(len=100) warnline
 
+  integer :: maxerrloc(3)
+  integer :: ii,jj,kk
+  CCTK_REAL :: maxerr
+
+
   CCTK_INT :: count
   CCTK_REAL :: s2, df, v2, w, vlowx, vlowy, vlowz 
   CCTK_REAL :: temp1, temp2, temp3
@@ -1401,8 +1410,8 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
   
   call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-    type2_bits = -1
-    excised = -1
+  type2_bits = -1
+  excised = -1
 
   nx = cctk_lsh(1)
   ny = cctk_lsh(2)
@@ -1440,18 +1449,27 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
         if ( (dens(i,j,k).le.sqrt(Whisky_Det(i,j,k)) * &
              whisky_rho_min*(1.0d0+whisky_atmo_tolerance)) &
              .or.(tau(i,j,k) .le. 0d0)) then
-          rho(i,j,k) = whisky_rho_min
-          dens(i,j,k) = sqrt(Whisky_det(i,j,k)) * rho(i,j,k)
-          scon(i,j,k,:) = 0.d0
-          w_lorentz(i,j,k) = 1.d0
-          press(i,j,k) = 0.1d0 * eosgeneral_pmin
-          eps(i,j,k) = press(i,j,k) / rho(i,j,k) ! Note that this should be improved
-          tau(i,j,k) = sqrt(Whisky_det(i,j,k)) * rho(i,j,k) * eps(i,j,k)
-!!$ call Whisky_DumpPointInAtmos(rho(i,j,k), &
-!!$ velx(i,j,k), vely(i,j,k), velz(i,j,k), & 
-!!$ eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & 
-!!$ dens(i,j,k), sx(i,j,k), sy(i,j,k), sz(i,j,k), tau(i,j,k), & 
-!!$  Whisky_det(i,j,k), whisky_rho_min)
+
+ !          write(6,*) "problem in atmo probably"
+ !          write(6,*) i,j,k
+ !          write(6,"(1P10E15.6)") dens(i,j,k), Whisky_Det(i,j,k), tau(i,j,k)
+ !          write(6,"(1P10E15.6)") eps(i,j,k)
+           rho(i,j,k) = whisky_rho_min
+           dens(i,j,k) = sqrt(Whisky_det(i,j,k)) * rho(i,j,k)
+           scon(i,j,k,:) = 0.d0
+           w_lorentz(i,j,k) = 1.d0
+           press(i,j,k) = 0.1d0 * eosgeneral_pmin
+!           eps(i,j,k) = press(i,j,k) / rho(i,j,k) ! Note that this should be improved
+           eps(i,j,k) = 1.0e-12
+           tau(i,j,k) = sqrt(Whisky_det(i,j,k)) * rho(i,j,k) * eps(i,j,k)
+
+           SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
+           atmosphere_mask(i,j,k) = 1
+!           call Whisky_DumpPointInAtmos(rho(i,j,k), &
+!                vel(i,j,k,1), vel(i,j,k,2), vel(i,j,k,3), & 
+!                eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & 
+!                dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), tau(i,j,k), & 
+!                Whisky_det(i,j,k), whisky_rho_min)
         end if
         
         udens = dens(i,j,k) / sqrt(Whisky_det(i,j,k))
@@ -1541,12 +1559,29 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
        if(whisky_reflevel.ge.whisky_c2p_warn_from_reflevel) then
           !$OMP CRITICAL
           call CCTK_WARN(1, 'Con2PrimGeneral: error: did not converge in ')
-          write(warnline,'(a20,g12.7,a10)') '              ',&
+          write(warnline,'(a20,i12,a10)') '              ',&
                whisky_countmax,' steps'
           call CCTK_WARN(1,warnline)
           write(warnline,'(a28,i2)') 'on carpet reflevel: ',whisky_reflevel
           call CCTK_WARN(1,warnline)
-          call CCTK_WARN(whisky_c2p_warnlevel, "Did not converge")
+          maxerr =  maxval(abs(press_new - press_old) & 
+               / abs(press_new))
+          maxerrloc = maxloc(abs(press_new - press_old) & 
+               / abs(press_new)) 
+          ii = maxerrloc(1)
+          jj = maxerrloc(2)
+          kk = maxerrloc(3)
+          write(warnline,'(1P10E15.6)') maxerr
+          call CCTK_WARN(1,warnline)
+          write(warnline,'(1P10E15.6)') press_new(ii,jj,kk), press_old(ii,jj,kk)
+          call CCTK_WARN(1,warnline)
+          write(warnline,'(3i4)') ii,jj,kk
+          call CCTK_WARN(1,warnline)
+          write(warnline,'(1P10E15.6)') x(ii,jj,kk), &
+               y(ii,jj,kk), &
+               z(ii,jj,kk)
+          call CCTK_WARN(1,warnline)
+          call CCTK_WARN(0,"Abort!")
           !$OMP END CRITICAL
           exit 
        else
@@ -1565,10 +1600,10 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
           !do not compute if in atmosphere
           !or in an excised region
 
-          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) then
-            press_new(i,j,k) = press(i,j,k)
-            press_old(i,j,k) = press(i,j,k)
-            cycle
+          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .ne. 0) &
+               then
+             press_new(i,j,k) = press(i,j,k)
+             press_old(i,j,k) = press(i,j,k)
           endif
 
 
@@ -1640,7 +1675,8 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
           point_loop_condition(i,j,k) = point_loop_condition(i,j,k) .and. &
                (abs(press_new(i,j,k) - press_old(i,j,k))  > &
                       whisky_del_ptol)          
-        end do
+          
+         end do
       end do
     end do
     
@@ -1654,7 +1690,7 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
     loop_condition = loop_condition .or. &
          (count < whisky_countmin) 
 
-  end do
+   end do
 
   deallocate(point_loop_condition)
 
@@ -1776,23 +1812,27 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
             call CCTK_WARN(1, '   specific internal energy just went below 0! ')
             write(warnline,'(a28,i2)') 'on carpet reflevel: ',whisky_reflevel
             call CCTK_WARN(1,warnline)
-            write(warnline,'(a20,3g16.7)') 'xyz location: ',&
+            write(warnline,'(a20,3i6)') 'ijk location: ',&
+                 i,j,k
+            call CCTK_WARN(1,warnline)
+            write(warnline,'(a20,1P10E15.6)') 'xyz location: ',&
                  x(i,j,k),y(i,j,k),z(i,j,k)
             call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,i4)') 'refinement level: ',whisky_reflevel
+            write(warnline,'(a25,i5)') 'refinement level: ',whisky_reflevel
             call CCTK_WARN(1,warnline)
-            write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k)
+            write(warnline,'(a20,1P1E15.6)') 'radius: ',r(i,j,k)
             call CCTK_WARN(1,warnline)
-            write(warnline,'(a20,3g16.7)') 'velocities: ',&
+            write(warnline,'(a20,1P1E15.6)') 'velocities: ',&
                  vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3)
             call CCTK_WARN(1,warnline)
-            call CCTK_WARN(whisky_c2p_warnlevel, \
+            call CCTK_WARN(whisky_c2p_warnlevel, &
             "Specific internal energy negative")
             !$OMP END CRITICAL
             exit
           else
             !$OMP CRITICAL
-            write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',whisky_reflevel
+            write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but \
+                        I was told to ignore level: ',whisky_reflevel
             call CCTK_WARN(1,warnline)
             write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
                  x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)