fixed surface tension bug

Author: sbryngelson

src/pre_process_code/m_global_parameters.f90

@@ -748,6 +748,8 @@ SUBROUTINE s_initialize_nonpoly
             pb0 = pb0/pl0
             pv = pv/pl0
 
+            print*, 'pb0 nondim/final', pb0
+
             ! bubble wall temperature, normalized by T0, in the liquid
             ! keeps a constant (cold liquid assumption)
             Tw = 1.D0

src/pre_process_code/m_initial_condition.f90

@@ -589,6 +589,7 @@ SUBROUTINE s_assign_patch_species_primitive_variables_bubbles( patch_id,j,k,l)
                 DO i = 1,nb
                     IF( q_prim_vf(bub_idx%ps(i))%sf(j,k,l) == dflt_real ) THEN
                         q_prim_vf(bub_idx%ps(i))%sf(j,k,l) = pb0(i) 
+                        print*, 'setting to pb0'
                     END IF
                     IF( q_prim_vf(bub_idx%ms(i))%sf(j,k,l) == dflt_real ) THEN
                         q_prim_vf(bub_idx%ms(i))%sf(j,k,l) = mass_v0(i)

src/simulation_code/m_bubbles.f90

@@ -183,7 +183,17 @@ SUBROUTINE s_compute_bubble_source(idir, q_prim_vf, q_cons_vf, mydivu, &
                            bub_m_src(q,j,k,l) = 0d0
                         END IF
                     END IF
+    
                 END DO; END DO; END DO; END DO
+
+                IF (DEBUG) THEN
+                    PRINT*, 'bub rhs'
+                    PRINT*, 'bub adv', bub_adv_src(:,0,0)
+                    PRINT*, 'bub r', bub_r_src(q,:,0,0)
+                    PRINT*, 'bub v', bub_v_src(q,:,0,0)
+                    PRINT*, 'bub p', bub_p_src(q,:,0,0)
+                    PRINT*, 'bub m', bub_m_src(q,:,0,0)
+                END IF
             END IF
 
         END SUBROUTINE s_compute_bubble_source
@@ -383,14 +393,17 @@ FUNCTION f_cpbw_KM( fR0, fR, fV, fpb )
             IF (polytropic) THEN
                 f_cpbw_KM = Ca*((fR0/fR)**(3.d0*gam)) - Ca + 1d0
                 IF (Web/=dflt_real) f_cpbw_KM = f_cpbw_KM + &
-                    (3.D0/(Web*fR0))*((fR0/fR)**(3.d0*gam))
+                    (2.D0/(Web*fR0))*((fR0/fR)**(3.d0*gam))
             ELSE
                 f_cpbw_KM = fpb 
+                ! @ t = 0, by default this is = pb0 = pl0[1] + 2*ss/(R0ref * R) computed by s_init_nonpoly 
             END IF
 
             ! PRINT*, 'surface tension component', (3.D0/(Web*fR0))*((fR0/fR)**(3.d0*gam))
 
-            IF (  Web /=dflt_real) f_cpbw_KM = f_cpbw_KM - 3.D0/(fR*Web)
+            IF (  Web /=dflt_real) f_cpbw_KM = f_cpbw_KM - 2.D0/(fR*Web)
+
+
             IF (Re_inv/=dflt_real) f_cpbw_KM = f_cpbw_KM - 4.D0*Re_inv*fV/fR
 
             ! PRINT*, ((fR0/fR)**(3.d0*gam))*(3.D0/(Web*fR0))-3.D0/(fR*Web)
@@ -435,6 +448,7 @@ FUNCTION f_rddot_KM( fpbdot, fCp, fCpbw, fRho, fR, fV, fR0, fC )
             tmp2 = 1.5D0*(fV**2d0)*( tmp1/2d0-1d0 ) +   &
                    (1d0 + tmp1)*(fCpbw - fCp)/fRho  +   &
                    cdot_star * fR/(fRho*fC)
+
 
             IF (Re_inv==dflt_real) THEN
                 f_rddot_KM = tmp2/( fR*(1d0-tmp1) )