Lagrange bubbles fixing PR part 5

Author: Diego Vaca

src/common/m_constants.fpp

@@ -46,25 +46,28 @@ module m_constants
     ! RKCK constants
     integer, parameter :: num_ts_rkck = 6 !< Number of time-stages in the RKCK stepper
     ! RKCK 4th/5th time stepper coefficients based on Cash J. and Karp A. (1990)
-    real(kind(0d0)), parameter :: rkck_c1 = 0.0d0, rkck_c2 = 0.2d0, rkck_c3 = 0.3d0, rkck_c4 = 0.6d0, &  ! c1 c2 c3 c4 c5 c6
+    real(kind(0d0)), parameter :: rkck_c1 = 0.0d0, rkck_c2 = 0.2d0, rkck_c3 = 0.3d0, rkck_c4 = 0.6d0, &     ! c1 c2 c3 c4 c5 c6
                                   rkck_c5 = 1.0d0, rkck_c6 = 0.875d0
-    real(kind(0.d0)), dimension(6), parameter :: rkck_coef1 = (/0.2d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0/)   ! a21
-    real(kind(0.d0)), dimension(6), parameter :: rkck_coef2 = (/3.0d0/40.0d0, 9.0d0/40.0d0, 0.0d0, 0.0d0, &  ! a31 a32
+    real(kind(0.d0)), dimension(6), parameter :: rkck_coef1 = (/0.2d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0/)  ! a21
+    real(kind(0.d0)), dimension(6), parameter :: rkck_coef2 = (/3.0d0/40.0d0, 9.0d0/40.0d0, 0.0d0, 0.0d0, & ! a31 a32
                                                                 0.0d0, 0.0d0/)
-    real(kind(0.d0)), dimension(6), parameter :: rkck_coef3 = (/0.3d0, -0.9d0, 1.2d0, 0.0d0, 0.0d0, 0.0d0/)  ! a41 a42 a43
-    real(kind(0.d0)), dimension(6), parameter :: rkck_coef4 = (/-11.0d0/54.0d0, 2.5d0, -70.0d0/27.0d0, &  ! a51 a52 a53 a54
+    real(kind(0.d0)), dimension(6), parameter :: rkck_coef3 = (/0.3d0, -0.9d0, 1.2d0, 0.0d0, 0.0d0, 0.0d0/) ! a41 a42 a43
+    real(kind(0.d0)), dimension(6), parameter :: rkck_coef4 = (/-11.0d0/54.0d0, 2.5d0, -70.0d0/27.0d0, &    ! a51 a52 a53 a54
                                                                 35.d0/27.d0, 0.0d0, 0.0d0/)
-    real(kind(0.d0)), dimension(6), parameter :: rkck_coef5 = (/1631.0d0/55296.0d0, 175.0d0/512.0d0, &  ! a61 a62 a63 a64 a65
-                                                                575.d0/13824.d0, 44275.d0/110592.d0, 253.d0/4096.d0, 0.0d0/)
-    real(kind(0.d0)), dimension(6), parameter :: rkck_coef6 = (/37.d0/378.d0, 0.0d0, 250.d0/621.d0, &  ! b1 b2 b3 b4 b5 b6
+    real(kind(0.d0)), dimension(6), parameter :: rkck_coef5 = (/1631.0d0/55296.0d0, 175.0d0/512.0d0, &      ! a61 a62 a63 a64 a65
+                                                                575.d0/13824.d0, 44275.d0/110592.d0, &
+                                                                253.d0/4096.d0, 0.0d0/)
+    real(kind(0.d0)), dimension(6), parameter :: rkck_coef6 = (/37.d0/378.d0, 0.0d0, 250.d0/621.d0, &       ! b1 b2 b3 b4 b5 b6
                                                                 125.0d0/594.0d0, 0.0d0, 512.0d0/1771.0d0/)
-    real(kind(0.d0)), dimension(6), parameter :: rkck_coefE = (/37.d0/378.d0 - 2825.0d0/27648.0d0, 0.0d0, &  ! er1 er2 er3 er4 er5 er6 (4th/5th error)
-                                                                250.d0/621.d0 - 18575.0d0/48384.0d0, 125.0d0/594.0d0 - 13525.0d0/55296.0d0, &
+    real(kind(0.d0)), dimension(6), parameter :: rkck_coefE = (/37.d0/378.d0 - 2825.0d0/27648.0d0, 0.0d0, & ! er1 er2 er3 er4 er5 er6 (4th/5th error)
+                                                                250.d0/621.d0 - 18575.0d0/48384.0d0, &
+                                                                125.0d0/594.0d0 - 13525.0d0/55296.0d0, &
                                                                 -277.0d0/14336.0d0, 512.0d0/1771.0d0 - 0.25d0/)
     ! Adaptive rkck constants
     real(kind(0d0)), parameter :: verysmall_dt = 1.d-14 !< Very small dt, stop stepper
     real(kind(0d0)), parameter :: SAFETY = 0.9d0 !< Next dt will be maximum 0.9*dt if truncation error is above tolerance.
-    real(kind(0d0)), parameter :: RNDDEC = 1.0d05 !< Need to round the relative truncation error (5th decimal) to avoid the inclusion of random decimals when dividing by a very small number (rkck_tolerance)
+    real(kind(0d0)), parameter :: RNDDEC = 1.0d05 !< Need to round the relative truncation error (5th decimal) to avoid the inclusion of random decimals
+    ! after dividing calculated tolerance by a very small number (rkck_tolerance)
     real(kind(0d0)), parameter :: PSHRNK = -0.25d0 !< Factor to reduce dt when truncation error above tolerance
     real(kind(0d0)), parameter :: SHRNKDT = 0.5d0 !< Factor to reduce dt due to negative bubble radius
     real(kind(0d0)), parameter :: ERRCON = 1.89d-4 !< Limit to slightly increase dt when truncation error is between ERRCON and 1

src/simulation/m_bubbles_EL.fpp

@@ -162,7 +162,7 @@ contains
 
     end subroutine s_initialize_bubbles_EL_module
 
-    !> The purpose of this procedure is to start lagrange bubble parameters applying nondimensionalization if needed
+    !> The purpose of this procedure is to start the lagrange bubble parameters applying nondimensionalization if needed
     subroutine s_start_lagrange_inputs()
 
         integer :: id_bubbles, id_host
@@ -197,7 +197,7 @@ contains
         Web = 1d0/ss
         Re_inv = mul0
 
-        ! Need improvements to accept polytropic gas compression, isothermal and adiabatic  thermal models, and
+        ! Need improvement to accept polytropic gas compression, isothermal and adiabatic thermal models, and
         ! the Gilmore and RP bubble models.
         polytropic = .false.    ! Forcing no polytropic model
         thermal = 3             ! Forcing constant transfer coefficient model based on Preston et al., 2007
@@ -350,7 +350,7 @@ contains
                                             rhol, gamma, pi_inf, qv, Re)
         dynP = 0d0
         do i = 1, num_dims
-            dynP = dynP + 0.5*q_cons_vf(contxe + i)%sf(cell(1), cell(2), cell(3))**2/rhol
+            dynP = dynP + 0.5d0*q_cons_vf(contxe + i)%sf(cell(1), cell(2), cell(3))**2/rhol
         end do
         pliq = (q_cons_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - dynP - pi_inf)/gamma
         if (pliq < 0) print *, "Negative pressure", proc_rank, &
@@ -359,14 +359,14 @@ contains
         ! Initial particle pressure
         gas_p(bub_id, 1) = pliq + 2d0*(1d0/Web)/bub_R0(bub_id)
         if ((1d0/Web) /= 0d0) then
-            pcrit = pv - 4d0*(1d0/Web)/(3.d0*sqrt(3d0*gas_p(bub_id, 1)*bub_R0(bub_id)**3/(2d0*(1d0/Web))))
+            pcrit = pv - 4d0*(1d0/Web)/(3d0*sqrt(3d0*gas_p(bub_id, 1)*bub_R0(bub_id)**3d0/(2d0*(1d0/Web))))
             pref = gas_p(bub_id, 1)
         else
             pcrit = 0d0
         end if
 
         ! Initial particle mass
-        volparticle = 4d0/3d0*pi*bub_R0(bub_id)**3 ! volume
+        volparticle = 4d0/3d0*pi*bub_R0(bub_id)**3d0 ! volume
         gas_mv(bub_id, 1) = pv*volparticle*(1d0/(R_v*Tw))*(massflag) ! vapermass
         gas_mg(bub_id) = (gas_p(bub_id, 1) - pv*(massflag))*volparticle*(1d0/(R_n*Tw)) ! gasmass
         if (gas_mg(bub_id) <= 0d0) then
@@ -473,9 +473,9 @@ contains
                 indomain = particle_in_domain(inputvals(1:3))
                 if (indomain .and. (id > 0)) then
                     bub_id = bub_id + 1
-                    nBubs = bub_id ! local number of bubbles
-                    lag_id(bub_id, 1) = id  !global ID
-                    lag_id(bub_id, 2) = bub_id  !local ID
+                    nBubs = bub_id                  ! local number of bubbles
+                    lag_id(bub_id, 1) = id          ! global ID
+                    lag_id(bub_id, 2) = bub_id      ! local ID
                     mtn_pos(bub_id, 1:3, 1) = inputvals(1:3)
                     mtn_posPrev(bub_id, 1:3, 1) = inputvals(4:6)
                     mtn_vel(bub_id, 1:3, 1) = inputvals(7:9)
@@ -540,14 +540,14 @@ contains
                 fV = intfc_vel(k, 2)
                 fpb = gas_p(k, 2)
                 pint = f_cpbw_KM(fR0, fR, fV, fpb)
-                pint = pint + 0.5d0*fV**2
+                pint = pint + 0.5d0*fV**2d0
                 if (lag_params%cluster_type == 2) then
                     bub_dphidt(k) = (paux - pint) + term2
                     ! Accounting for the potential induced by the bubble averaged over the control volume
                     ! Note that this is based on the incompressible flow assumption near the bubble.
                     Rb = intfc_rad(k, 2)
                     term1_fac = 3d0/2d0*(Rb*(Romega**2d0 - Rb**2d0))/(Romega**3d0 - Rb**3d0)
-                    bub_dphidt(k) = bub_dphidt(k)/(1 - term1_fac)
+                    bub_dphidt(k) = bub_dphidt(k)/(1d0 - term1_fac)
                 end if
             end do
         end if
@@ -594,10 +594,10 @@ contains
                 end do
                 call s_convert_species_to_mixture_variables_acc(rhol, gamma, pi_inf, qv, myalpha, &
                                                                 myalpha_rho, Re, cell(1), cell(2), cell(3))
-                call s_compute_cson_from_pinf(k, q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
+                call s_compute_cson_from_pinf(q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
 
                 ! Velocity correction due to massflux
-                velint = fV - gas_dmvdt(k, stage)/(4d0*pi*fR**2*rhol)
+                velint = fV - gas_dmvdt(k, stage)/(4d0*pi*fR**2d0*rhol)
 
                 ! Interphase acceleration and update vars
                 intfc_dveldt(k, stage) = f_rddot_KM(fpbdt, pinf, pliqint, rhol, fR, velint, fR0, cson)
@@ -641,20 +641,18 @@ contains
         !! @param gamma Liquid specific heat ratio
         !! @param pi_inf Liquid stiffness
         !! @param cson Calculated speed of sound
-    subroutine s_compute_cson_from_pinf(bub_id, q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
+    subroutine s_compute_cson_from_pinf(q_prim_vf, pinf, cell, rhol, gamma, pi_inf, cson)
 #ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_cson_from_pinf
 #else
         !$acc routine seq
 #endif
-        integer, intent(in) :: bub_id
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
         real(kind(0d0)), intent(in) :: pinf, rhol, gamma, pi_inf
         integer, dimension(3), intent(in) :: cell
         real(kind(0d0)), intent(out) :: cson
 
         real(kind(0d0)) :: E, H
-        real(kind(0d0)), dimension(3) :: scoord
         real(kind(0d0)), dimension(num_dims) :: vel
         integer :: i
 
@@ -817,7 +815,7 @@ contains
         real(kind(0d0)), intent(out), optional :: preterm1, term2, Romega
 
         real(kind(0d0)), dimension(3) :: scoord, psi
-        real(kind(0d0)) :: dc, vol, aux, dist
+        real(kind(0d0)) :: dc, vol, aux
         real(kind(0d0)) :: volgas, term1, Rbeq, denom
         real(kind(0d0)) :: charvol, charpres, charvol2, charpres2
         integer, dimension(3) :: cellaux
@@ -998,22 +996,22 @@ contains
 
         end if
 
-        if (lag_params%pressure_corrector .and. present(preterm1)) then
+        if (lag_params%pressure_corrector) then
 
             !Valid if only one bubble exists per cell
-            volgas = intfc_rad(bub_id, 2)**3
-            denom = intfc_rad(bub_id, 2)**2
-            term1 = bub_dphidt(bub_id)*intfc_rad(bub_id, 2)**2
-            term2 = intfc_vel(bub_id, 2)*intfc_rad(bub_id, 2)**2
+            volgas = intfc_rad(bub_id, 2)**3d0
+            denom = intfc_rad(bub_id, 2)**2d0
+            term1 = bub_dphidt(bub_id)*intfc_rad(bub_id, 2)**2d0
+            term2 = intfc_vel(bub_id, 2)*intfc_rad(bub_id, 2)**2d0
 
             Rbeq = volgas**(1d0/3d0) !surrogate bubble radius
-            aux = dc**3 - Rbeq**3
+            aux = dc**3d0 - Rbeq**3d0
             term2 = term2/denom
-            term2 = 3d0/2d0*term2**2*Rbeq**3*(1d0 - Rbeq/dc)/aux
-            preterm1 = 3d0/2d0*Rbeq*(dc**2 - Rbeq**2)/(aux*denom)
+            term2 = 3d0/2d0*term2**2d0*Rbeq**3d0*(1d0 - Rbeq/dc)/aux
+            preterm1 = 3d0/2d0*Rbeq*(dc**2d0 - Rbeq**2d0)/(aux*denom)
 
             !Control volume radius
-            if (present(Romega)) Romega = dc
+            if (ptype == 2) Romega = dc
 
             ! Getting p_inf
             if (ptype == 1) then
@@ -1167,7 +1165,7 @@ contains
 
         lag_largestep = 0d0
         remove_id = 0
-        !$acc parallel loop gang vector default(present) reduction(+: lag_largestep, remove_id) private(k) copyin(RKstep)
+        !$acc parallel loop gang vector default(present) reduction(+: lag_largestep) reduction(MAX: remove_id) private(k) copyin(RKstep)
         do k = 1, nBubs
 
             radiusOld = intfc_rad(k, 2)
@@ -1188,7 +1186,7 @@ contains
                 print *, 'Negative bubble radius encountered'
                 lag_largestep = lag_largestep + 1d0
                 if (dt < 2d0*verysmall_dt) then
-                    remove_id = k
+                    remove_id = max(remove_id, k)
                 end if
             end if
 
@@ -1382,9 +1380,9 @@ contains
           !! @param scoord Calculated particle coordinates
     subroutine s_locate_cell(pos, cell, scoord)
 
-        real(kind(0d0)), dimension(3) :: pos
-        real(kind(0d0)), dimension(3), optional :: scoord
-        integer, dimension(3) :: cell
+        real(kind(0d0)), dimension(3), intent(in) :: pos
+        real(kind(0d0)), dimension(3), intent(out) :: scoord
+        integer, dimension(3), intent(inout) :: cell
 
         integer :: i
 
@@ -1420,16 +1418,14 @@ contains
         ! In other words,  the coordinate of the center of the cell is x_cc(cell).
 
         !coordinates in computational space
-        if (present(scoord)) then
-            scoord(1) = cell(1) + (pos(1) - x_cb(cell(1) - 1))/dx(cell(1))
-            scoord(2) = cell(2) + (pos(2) - y_cb(cell(2) - 1))/dy(cell(2))
-            scoord(3) = 0d0
-            if (p > 0) scoord(3) = cell(3) + (pos(3) - z_cb(cell(3) - 1))/dz(cell(3))
-            cell(:) = int(scoord(:))
-            do i = 1, num_dims
-                if (scoord(i) < 0.0d0) cell(i) = cell(i) - 1
-            end do
-        end if
+        scoord(1) = cell(1) + (pos(1) - x_cb(cell(1) - 1))/dx(cell(1))
+        scoord(2) = cell(2) + (pos(2) - y_cb(cell(2) - 1))/dy(cell(2))
+        scoord(3) = 0d0
+        if (p > 0) scoord(3) = cell(3) + (pos(3) - z_cb(cell(3) - 1))/dz(cell(3))
+        cell(:) = int(scoord(:))
+        do i = 1, num_dims
+            if (scoord(i) < 0.0d0) cell(i) = cell(i) - 1
+        end do
 
     end subroutine s_locate_cell
 
@@ -1847,7 +1843,7 @@ contains
 
         integer :: k
 
-        !$acc parallel loop gang vector default(present) reduction(MAX: Rmax_glb) reduction(MAX: Rmin_glb) private(k)
+        !$acc parallel loop gang vector default(present) reduction(MAX: Rmax_glb) reduction(MIN: Rmin_glb) private(k)
         do k = 1, nBubs
             Rmax_glb = max(Rmax_glb, intfc_rad(k, 1)/bub_R0(k))
             Rmin_glb = min(Rmin_glb, intfc_rad(k, 1)/bub_R0(k))

src/simulation/m_bubbles_EL_kernels.fpp

@@ -62,12 +62,12 @@ contains
         !$acc parallel loop gang vector default(present) private(l, s_coord, cell)
         do l = 1, nBubs
 
-            volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3
+            volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3d0
             s_coord(1:3) = lbk_s(l, 1:3, 2)
             call s_get_cell(s_coord, cell)
 
             strength_vol = volpart
-            strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2*lbk_vel(l, 2)
+            strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2d0*lbk_vel(l, 2)
 
             if (num_dims == 2) then
                 Vol = dx(cell(1))*dy(cell(2))*lag_params%charwidth
@@ -128,15 +128,15 @@ contains
         do l = 1, nBubs
             nodecoord(1:3) = 0
             center(1:3) = 0.0d0
-            volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3
+            volpart = 4.0d0/3.0d0*pi*lbk_rad(l, 2)**3d0
             s_coord(1:3) = lbk_s(l, 1:3, 2)
             center(1:2) = lbk_pos(l, 1:2, 2)
             if (p > 0) center(3) = lbk_pos(l, 3, 2)
             call s_get_cell(s_coord, cell)
             call s_compute_stddsv(cell, volpart, stddsv)
 
             strength_vol = volpart
-            strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2*lbk_vel(l, 2)
+            strength_vel = 4.0d0*pi*lbk_rad(l, 2)**2d0*lbk_vel(l, 2)
 
             !$acc loop collapse(3) private(cellaux, nodecoord)
             do i = 1, smearGrid
@@ -218,25 +218,25 @@ contains
         real(kind(0.d0)), intent(out) :: func
 
         real(kind(0.d0)) :: distance
-        real(kind(0.d0)) :: theta, dtheta, L2, dz, Lz2
+        real(kind(0.d0)) :: theta, dtheta, L2, dzp, Lz2
         integer :: Nr, Nr_count
 
-        distance = sqrt((center(1) - nodecoord(1))**2 + (center(2) - nodecoord(2))**2 + (center(3) - nodecoord(3))**2)
+        distance = sqrt((center(1) - nodecoord(1))**2d0 + (center(2) - nodecoord(2))**2d0 + (center(3) - nodecoord(3))**2d0)
 
         if (num_dims == 3) then
             !< 3D gaussian function
-            func = exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**3
+            func = exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**3d0
         else
             if (cyl_coord) then
                 !< 2D cylindrical function:
                 ! We smear particles in the azimuthal direction for given r
                 theta = 0d0
-                Nr = ceiling(2d0*PI*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
-                dtheta = 2d0*PI/Nr
+                Nr = ceiling(2d0*pi*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
+                dtheta = 2d0*pi/Nr
                 L2 = center(2)**2d0 + nodecoord(2)**2d0 - 2d0*center(2)*nodecoord(2)*cos(theta)
                 distance = DSQRT((center(1) - nodecoord(1))**2d0 + L2)
                 ! Factor 2d0 is for symmetry (upper half of the 2D field (+r) is considered)
-                func = dtheta/2d0/PI*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**3
+                func = dtheta/2d0/pi*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**3d0
                 Nr_count = 0
                 do while (Nr_count < Nr - 1)
                     Nr_count = Nr_count + 1
@@ -246,7 +246,7 @@ contains
                     distance = DSQRT((center(1) - nodecoord(1))**2d0 + L2)
                     ! nodecoord(2)*dtheta is the azimuthal width of the cell
                     func = func + &
-                           dtheta/2d0/PI*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**(3*(strength_idx + 1))
+                           dtheta/2d0/pi*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**(3d0*(strength_idx + 1))
                 end do
             else
 
@@ -255,16 +255,16 @@ contains
                 theta = 0d0
                 Nr = ceiling(lag_params%charwidth/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
                 Nr_count = 1 - mapCells
-                dz = y_cb(cellaux(2) + 1) - y_cb(cellaux(2))
-                Lz2 = (center(3) - (dz*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
+                dzp = y_cb(cellaux(2) + 1) - y_cb(cellaux(2))
+                Lz2 = (center(3) - (dzp*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
                 distance = DSQRT((center(1) - nodecoord(1))**2d0 + (center(2) - nodecoord(2))**2d0 + Lz2)
-                func = dz/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**3
+                func = dzp/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**3d0
                 do while (Nr_count < Nr - 1 + (mapCells - 1))
                     Nr_count = Nr_count + 1
-                    Lz2 = (center(3) - (dz*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
+                    Lz2 = (center(3) - (dzp*(5d-1 + Nr_count) - lag_params%charwidth/2d0))**2d0
                     distance = DSQRT((center(1) - nodecoord(1))**2d0 + (center(2) - nodecoord(2))**2d0 + Lz2)
                     func = func + &
-                           dz/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2)/(DSQRT(2.0d0*pi)*stddsv)**(3*(strength_idx + 1))
+                           dzp/lag_params%charwidth*exp(-0.5d0*(distance/stddsv)**2d0)/(DSQRT(2.0d0*pi)*stddsv)**(3d0*(strength_idx + 1))
                 end do
             end if
         end if

src/simulation/m_global_parameters.fpp

@@ -699,7 +699,7 @@ contains
         lag_params%write_bubbles = .false.
         lag_params%write_bubbles_stats = .false.
         lag_params%nBubs_glb = dflt_int
-        lag_params%epsilonb = dflt_real
+        lag_params%epsilonb = 1.0d0
         lag_params%charwidth = dflt_real
         lag_params%valmaxvoid = dflt_real
         lag_params%c0 = dflt_real