minor changes

Author: sbryngelson

src/simulation_code/m_global_parameters.f90

@@ -1148,26 +1148,23 @@ SUBROUTINE s_comp_n_from_cons( vftmp,nRtmp,ntmp )
             CALL s_quad( nRtmp**3d0,nR3 )
 
             IF ( nR3 < 0d0 ) THEN
-                PRINT*, vftmp, nR3, nRtmp(:)
                 ! DO i = 1,nb
                     ! IF (nRtmp(i) < small_alf) THEN
                         ! nRtmp(i) = small_alf
                     ! END IF
                 ! END DO
-                STOP 'nR3 is negative'
-                nR3 = 1.d-12
+                ! nR3 = 1.d-12
+                PRINT*, vftmp, nR3, nRtmp(:)
+                STOP'nR3 is negative'
             END IF
             IF (vftmp < 0d0) THEN
-                PRINT*, vftmp, nR3, nRtmp(:)
                 ! vftmp = small_alf
+                ! ntmp = DSQRT( (4.d0*pi/3.d0)*nR3/1.d-12 )
+                PRINT*, vftmp, nR3, nRtmp(:)
                 STOP 'vf negative'
-                ntmp = DSQRT( (4.d0*pi/3.d0)*nR3/1.d-12 )
-            ELSE
-                ntmp = DSQRT( (4.d0*pi/3.d0)*nR3/vftmp )
             END IF
 
-
-            ! ntmp = 1d0
+            ntmp = DSQRT( (4.d0*pi/3.d0)*nR3/vftmp )
 
         END SUBROUTINE s_comp_n_from_cons
 
@@ -1184,6 +1181,16 @@ SUBROUTINE s_comp_n_from_prim( vftmp,Rtmp,ntmp )
             REAL(KIND(0.D0)) :: R3
 
             CALL s_quad( Rtmp**3d0,R3 ) 
+
+            IF ( R3 < 0d0 ) THEN
+                PRINT*, vftmp, R3, Rtmp(:)
+                STOP 'R3 is negative'
+            END IF
+            IF (vftmp < 0d0) THEN
+                PRINT*, vftmp, R3, Rtmp(:)
+                STOP 'vf negative'
+            END IF
+
             ntmp = (3.d0/(4.d0*pi)) * vftmp/R3
             ! ntmp = 1d0
 

src/simulation_code/m_qbmm.f90

@@ -178,8 +178,8 @@ SUBROUTINE s_coeff( pres,rho,c,coeffs )
                         coeffs(2,i1,i2) = -3d0*i2/2d0
                         coeffs(3,i1,i2) = i2/rho
                         coeffs(4,i1,i2) = i1
-                        IF (Re_inv/=dflt_real) coeffs(5,i1,i2) = 4d0*i2*Re_inv/rho
-                        IF (  Web /=dflt_real) coeffs(6,i1,i2) = 2*i2/Web/rho
+                        IF (Re_inv/=dflt_real) coeffs(5,i1,i2) = -4d0*i2*Re_inv/rho
+                        IF (  Web /=dflt_real) coeffs(6,i1,i2) = -2d0*i2/Web/rho
                     ELSE IF (bubble_model==2) THEN
                         ! KM with approximation of 1/(1-V/C) = 1+V/C
                         coeffs(1,i1,i2)  = -3d0*i2/2d0
@@ -228,7 +228,12 @@ SUBROUTINE s_mom_inv( q_prim_vf, momsp, moms3d, is1, is2, is3 )
                     n_tait = fluid_pp(1)%gamma
                     n_tait = 1.d0/n_tait + 1.d0 !make this the usual little 'gamma'
                     B_tait = fluid_pp(1)%pi_inf
-                    c =  DSQRT(n_tait*(pres+B_tait)/(rho*(1.d0-alf)))
+                    c = n_tait*(pres+B_tait)/(rho*(1.d0-alf))
+                    IF (c > 0.d0) THEN
+                        c = DSQRT(c)
+                    ELSE
+                        c = sgm_eps
+                    END IF                    
                 END IF
 
                 CALL s_coeff(pres,rho,c,coeff)
@@ -284,15 +289,25 @@ SUBROUTINE s_mom_inv( q_prim_vf, momsp, moms3d, is1, is2, is3 )
                     momsp(1)%sf(id1,id2,id3) = f_quad(abscX,abscY,wght,3d0,0d0,0d0)
                     momsp(2)%sf(id1,id2,id3) = 4.d0*pi*nbub*f_quad(abscX,abscY,wght,2d0,1d0,0d0)
                     momsp(3)%sf(id1,id2,id3) = f_quad(abscX,abscY,wght,3d0,2d0,0d0)
-                    momsp(4)%sf(id1,id2,id3) = f_quad(abscX,abscY,wght,3d0*(1d0-gam),0d0,3d0*gam)
-
-                    ! DO i1 = 1,nterms
-                    !     IF (momsp(i1)%sf(id1,id2,id3)/=momsp(i1)%sf(id1,id2,id3)) THEN
-                    !         PRINT*, 'nan in momsp', i1,id1
-                    !         PRINT*, 'moms: ', moms(:)
-                    !         CALL s_mpi_abort()
-                    !     END IF
-                    ! END DO
+                    IF (ABS(gam-1.d0) <= 1.d-4) THEN
+                        ! Gam \approx 1, don't risk imaginary quadrature
+                        momsp(4)%sf(id1,id2,id3) = 1.d0
+                    ELSE 
+                        momsp(4)%sf(id1,id2,id3) = f_quad(abscX,abscY,wght,3d0*(1d0-gam),0d0,3d0*gam)
+                    END IF
+
+                    DO i1 = 1,4
+                        IF (momsp(i1)%sf(id1,id2,id3) /= momsp(i1)%sf(id1,id2,id3)) THEN
+                            PRINT*, 'NaN in sp moment', i1, 'location',id1,id2,id3
+                            PRINT*, 'Rs', Rvec(:)
+                            PRINT*, 'alpha', alf
+                            PRINT*, 'nbub', nbub
+                            PRINT*, 'abscX', abscX(:,:)
+                            PRINT*, 'abscY', abscY(:,:)
+                            PRINT*, 'wght' , wght(:,:)
+                            CALL s_mpi_abort()
+                        END IF
+                    END DO
                 ELSE
                     DO q = 1,nb
                         DO i1 = 0,2; DO i2 = 0,2
@@ -384,7 +399,7 @@ SUBROUTINE s_hyqmom(frho,fup,fmom)
             c2 = d2 - bu**2d0 
             frho(1) = fmom(1)/2d0; 
             frho(2) = fmom(1)/2d0; 
-            c2 = MAXVAL( (/ c2,verysmall /) )
+            c2 = MAXVAL( (/ c2, verysmall /) )
             fup(1) = bu - DSQRT(c2)
             fup(2) = bu + DSQRT(c2) 
 

src/simulation_code/m_riemann_solvers.f90

@@ -2743,9 +2743,9 @@ SUBROUTINE s_compute_arithmetic_average_state(j,k,l) ! --------
                 END IF
 
                 IF ( (ptilde_L .ne. ptilde_L) .or. (ptilde_R .ne. ptilde_R) ) THEN
-                    PRINT*, 'Ptilde NaN'
-                    PRINT*, alpha_L, pres_L, PbwR3Lbar, R3Lbar, rho_L, R3V2Lbar, R3Lbar
-                    PRINT*, alpha_R, pres_R, PbwR3Rbar, R3Rbar, rho_R, R3V2Rbar, R3Rbar
+                    PRINT*, 'Ptilde NaN at ',j,k,l,x_cb(j)
+                    PRINT*, nbub_L, alpha_L, pres_L, PbwR3Lbar, R3Lbar, rho_L, R3V2Lbar, R3Lbar
+                    PRINT*, nbub_R, alpha_R, pres_R, PbwR3Rbar, R3Rbar, rho_R, R3V2Rbar, R3Rbar
                     CALL s_mpi_abort()
                 END IF