Add files necessary for compile from MovingBubblesFresh

This commit adds files necessary for compiling and running the EL particle solver. These files/changes are not part of this specific feature development. These changes are part of the MovingBubblesFresh commits.

Since this feature for moving solid lagrangian particles is a fork of the MovingBubblesFresh development fork, it has inherent dependencies. When MovingBubblesFresh gets merged into MFC, this commit will be irrelevant.

Author: jaguilar

src/simulation/m_bubbles.fpp

@@ -17,6 +17,8 @@ module m_bubbles
 
     use m_helper_basic         !< Functions to compare floating point numbers
 
+    use m_bubbles_EL_kernels
+
     implicit none
 
     real(wp) :: chi_vw  !< Bubble wall properties (Ando 2010)
@@ -40,7 +42,7 @@ contains
         !!  @param f_bub_adv_src Source for bubble volume fraction
         !!  @param f_divu Divergence of velocity
         !!  @param fCson Speed of sound from fP (EL)
-    elemental function f_rddot(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, fntait, fBtait, f_bub_adv_src, f_divu, fCson)
+    function f_rddot(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, fntait, fBtait, f_bub_adv_src, f_divu, fCson)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
         real(wp), intent(in) :: fntait, fBtait, f_bub_adv_src, f_divu
@@ -84,7 +86,7 @@ contains
         !!  @param fR Current bubble radius
         !!  @param fV Current bubble velocity
         !!  @param fpb Internal bubble pressure
-    elemental function f_cpbw(fR0, fR, fV, fpb)
+    function f_cpbw(fR0, fR, fV, fpb)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fR0, fR, fV, fpb
 
@@ -103,7 +105,7 @@ contains
         !!  @param fCpinf Driving bubble pressure
         !!  @param fntait Tait EOS parameter
         !!  @param fBtait Tait EOS parameter
-    elemental function f_H(fCpbw, fCpinf, fntait, fBtait)
+    function f_H(fCpbw, fCpinf, fntait, fBtait)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fCpbw, fCpinf, fntait, fBtait
 
@@ -123,7 +125,7 @@ contains
         !! @param fntait Tait EOS parameter
         !! @param fBtait Tait EOS parameter
         !! @param fH Bubble enthalpy
-    elemental function f_cgas(fCpinf, fntait, fBtait, fH)
+    function f_cgas(fCpinf, fntait, fBtait, fH)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fCpinf, fntait, fBtait, fH
 
@@ -146,7 +148,7 @@ contains
         !!  @param fBtait Tait EOS parameter
         !!  @param advsrc Advection equation source term
         !!  @param divu Divergence of velocity
-    elemental function f_cpinfdot(fRho, fP, falf, fntait, fBtait, advsrc, divu)
+    function f_cpinfdot(fRho, fP, falf, fntait, fBtait, advsrc, divu)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fRho, fP, falf, fntait, fBtait, advsrc, divu
 
@@ -176,7 +178,7 @@ contains
         !!  @param fV Current bubble velocity
         !!  @param fR0 Equilibrium bubble radius
         !!  @param fpbdot Time derivative of the internal bubble pressure
-    elemental function f_Hdot(fCpbw, fCpinf, fCpinf_dot, fntait, fBtait, fR, fV, fR0, fpbdot)
+    function f_Hdot(fCpbw, fCpinf, fCpinf_dot, fntait, fBtait, fR, fV, fR0, fpbdot)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fCpbw, fCpinf, fCpinf_dot, fntait, fBtait
         real(wp), intent(in) :: fR, fV, fR0, fpbdot
@@ -211,7 +213,7 @@ contains
         !!  @param fR Current bubble radius
         !!  @param fV Current bubble velocity
         !!  @param fCpbw Boundary wall pressure
-    elemental function f_rddot_RP(fCp, fRho, fR, fV, fCpbw)
+    function f_rddot_RP(fCp, fRho, fR, fV, fCpbw)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fCp, fRho, fR, fV, fCpbw
 
@@ -234,7 +236,7 @@ contains
         !!  @param fcgas Current gas sound speed
         !!  @param fntait Tait EOS parameter
         !!  @param fBtait Tait EOS parameter
-    elemental function f_rddot_G(fCpbw, fR, fV, fH, fHdot, fcgas, fntait, fBtait)
+    function f_rddot_G(fCpbw, fR, fV, fH, fHdot, fcgas, fntait, fBtait)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fCpbw, fR, fV, fH, fHdot
         real(wp), intent(in) :: fcgas, fntait, fBtait
@@ -257,7 +259,7 @@ contains
         !!  @param fR Current bubble radius
         !!  @param fV Current bubble velocity
         !!  @param fpb Internal bubble pressure
-    elemental function f_cpbw_KM(fR0, fR, fV, fpb)
+    function f_cpbw_KM(fR0, fR, fV, fpb)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fR0, fR, fV, fpb
         real(wp) :: f_cpbw_KM
@@ -284,7 +286,7 @@ contains
         !!  @param fV Current bubble velocity
         !!  @param fR0 Equilibrium bubble radius
         !!  @param fC Current sound speed
-    elemental function f_rddot_KM(fpbdot, fCp, fCpbw, fRho, fR, fV, fR0, fC)
+    function f_rddot_KM(fpbdot, fCp, fCpbw, fRho, fR, fV, fR0, fC)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fpbdot, fCp, fCpbw
         real(wp), intent(in) :: fRho, fR, fV, fR0, fC
@@ -318,7 +320,7 @@ contains
     !>  Subroutine that computes bubble wall properties for vapor bubbles
         !!  @param pb_in Internal bubble pressure
         !!  @param iR0 Current bubble size index
-    elemental subroutine s_bwproperty(pb_in, iR0, chi_vw_out, k_mw_out, rho_mw_out)
+    subroutine s_bwproperty(pb_in, iR0, chi_vw_out, k_mw_out, rho_mw_out)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: pb_in
         integer, intent(in) :: iR0
@@ -349,7 +351,7 @@ contains
         !!  @param fbeta_c Mass transfer coefficient (EL)
         !!  @param fR_m Mixture gas constant (EL)
         !!  @param fgamma_m Mixture gamma (EL)
-    elemental subroutine s_vflux(fR, fV, fpb, fmass_v, iR0, vflux, fmass_g, fbeta_c, fR_m, fgamma_m)
+    subroutine s_vflux(fR, fV, fpb, fmass_v, iR0, vflux, fmass_g, fbeta_c, fR_m, fgamma_m)
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fR
         real(wp), intent(in) :: fV
@@ -407,7 +409,8 @@ contains
         !!  @param fbeta_t Mass transfer coefficient (EL)
         !!  @param fR_m Mixture gas constant (EL)
         !!  @param fgamma_m Mixture gamma (EL)
-    elemental function f_bpres_dot(fvflux, fR, fV, fpb, fmass_v, iR0, fbeta_t, fR_m, fgamma_m)
+    function f_bpres_dot(fvflux, fR, fV, fpb, fmass_v, iR0, fbeta_t, fR_m, fgamma_m)
+        !$DIR INLINENEVER f_bpres_dot
         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in) :: fvflux
         real(wp), intent(in) :: fR
@@ -463,27 +466,32 @@ contains
         !!  @param fbeta_t Heat transfer coefficient (EL)
         !!  @param fCson Speed of sound (EL)
         !!  @param adap_dt_stop Fail-safe exit if max iteration count reached
-    subroutine s_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
-                              fntait, fBtait, f_bub_adv_src, f_divu, &
-                              bub_id, fmass_v, fmass_g, fbeta_c, &
-                              fbeta_t, fCson, adap_dt_stop)
-        $:GPU_ROUTINE(function_name='s_advance_step',parallelism='[seq]', &
-            & cray_inline=True)
+    function f_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+                            fntait, fBtait, f_bub_adv_src, f_divu, &
+                            bub_id, fmass_v, fmass_g, fbeta_c, &
+                            fbeta_t, fCson, fRe, fPos, &
+                            fVel, cell, q_prim_vf) result(adap_dt_stop)
+        $:GPU_ROUTINE(parallelism='[seq]')
 
         real(wp), intent(inout) :: fR, fV, fpb, fmass_v
         real(wp), intent(in) :: fRho, fP, fR0, fpbdot, alf
         real(wp), intent(in) :: fntait, fBtait, f_bub_adv_src, f_divu
         integer, intent(in) :: bub_id
         real(wp), intent(in) :: fmass_g, fbeta_c, fbeta_t, fCson
-        integer, intent(inout) :: adap_dt_stop
+        real(wp), intent(inout), dimension(3), optional :: fPos, fVel
+        real(wp), intent(in), optional :: fRe
+        integer, intent(in), dimension(3), optional :: cell
+        type(scalar_field), intent(in), dimension(sys_size), optional :: q_prim_vf
 
         real(wp), dimension(5) :: err !< Error estimates for adaptive time stepping
         real(wp) :: t_new !< Updated time step size
         real(wp) :: h0, h !< Time step size
         real(wp), dimension(4) :: myR_tmp1, myV_tmp1, myR_tmp2, myV_tmp2 !< Bubble radius, radial velocity, and radial acceleration for the inner loop
         real(wp), dimension(4) :: myPb_tmp1, myMv_tmp1, myPb_tmp2, myMv_tmp2 !< Gas pressure and vapor mass for the inner loop (EL)
-        real(wp) :: fR2, fV2, fpb2, fmass_v2
-        integer :: iter_count
+        real(wp) :: fR2, fV2, fpb2, fmass_v2, f_bTemp
+        real(wp), dimension(3) :: vTemp, aTemp
+        integer :: adap_dt_stop
+        integer :: l, iter_count
 
         call s_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
                                  fntait, fBtait, f_bub_adv_src, f_divu, fCson, h0)
@@ -565,6 +573,42 @@ contains
                         ! Update pb and mass_v
                         fpb = myPb_tmp1(4)
                         fmass_v = myMv_tmp1(4)
+
+                        select case (lag_vel_model)
+                        case (1)
+                            do l = 1, num_dims
+                                vTemp(l) = f_interpolate_velocity(fR, cell, l, q_prim_vf)
+                            end do
+                            do l = 1, num_dims
+                                fVel(l) = vTemp(l)
+                                fPos(l) = fPos(l) + h*vTemp(l)
+                            end do
+                        case (2)
+                            do l = 1, num_dims
+                                f_bTemp = f_get_bubble_force(fPos(l), fR, fV, fVel(l), fmass_g, fmass_v, &
+                                                             fRe, fRho, cell, l, q_prim_vf)
+                                aTemp(l) = f_bTemp/(fmass_g + fmass_v)
+                            end do
+                            do l = 1, num_dims
+                                fVel(l) = fVel(l) + h*aTemp(l)
+                                fPos(l) = fPos(l) + h*fVel(l)
+                            end do
+                        case (3)
+                            do l = 1, num_dims
+                                f_bTemp = f_get_bubble_force(fPos(l), fR, fV, fVel(l), fmass_g, fmass_v, &
+                                                             fRe, fRho, cell, l, q_prim_vf)
+                                aTemp(l) = 2._wp*f_bTemp/(fmass_g + fmass_v) - 3*fV*fVel(l)/fR
+                            end do
+                            do l = 1, num_dims
+                                fVel(l) = fVel(l) + h*aTemp(l)
+                                fPos(l) = fPos(l) + h*fVel(l)
+                            end do
+                        case default
+                            do l = 1, num_dims
+                                fVel(l) = fVel(l)
+                                fPos(l) = fPos(l)
+                            end do
+                        end select
                     end if
 
                     ! Update step size for the next sub-step
@@ -588,7 +632,7 @@ contains
 
         if (iter_count >= adap_dt_max_iters) adap_dt_stop = 1
 
-    end subroutine s_advance_step
+    end function f_advance_step
 
     !> Choose the initial time step size for the adaptive time stepping routine
         !!  (See Heirer, E. Hairer S.P.Nørsett G. Wanner, Solving Ordinary
@@ -613,10 +657,10 @@ contains
         $:GPU_ROUTINE(function_name='s_initial_substep_h',parallelism='[seq]', &
             & cray_inline=True)
 
-        real(wp), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
-        real(wp), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu
-        real(wp), intent(IN) :: fCson
-        real(wp), intent(OUT) :: h
+        real(wp), intent(in) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
+        real(wp), intent(in) :: fntait, fBtait, f_bub_adv_src, f_divu
+        real(wp), intent(in) :: fCson
+        real(wp), intent(out) :: h
 
         real(wp), dimension(2) :: h_size !< Time step size (h0, h1)
         real(wp), dimension(3) :: d_norms !< norms (d_0, d_1, d_2)
@@ -697,12 +741,12 @@ contains
         $:GPU_ROUTINE(function_name='s_advance_substep',parallelism='[seq]', &
             & cray_inline=True)
 
-        real(wp), intent(OUT) :: err
-        real(wp), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
-        real(wp), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu, h
-        integer, intent(IN) :: bub_id
-        real(wp), intent(IN) :: fmass_v, fmass_g, fbeta_c, fbeta_t, fCson
-        real(wp), dimension(4), intent(OUT) :: myR_tmp, myV_tmp, myPb_tmp, myMv_tmp
+        real(wp), intent(out) :: err
+        real(wp), intent(in) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
+        real(wp), intent(in) :: fntait, fBtait, f_bub_adv_src, f_divu, h
+        integer, intent(in) :: bub_id
+        real(wp), intent(in) :: fmass_v, fmass_g, fbeta_c, fbeta_t, fCson
+        real(wp), dimension(4), intent(out) :: myR_tmp, myV_tmp, myPb_tmp, myMv_tmp
 
         real(wp), dimension(4) :: myA_tmp, mydPbdt_tmp, mydMvdt_tmp
         real(wp) :: err_R, err_V
@@ -802,14 +846,15 @@ contains
         !!  @param fbeta_c Mass transfer coefficient
         !!  @param fbeta_t Heat transfer coefficient
         !!  @param fdPbdt_tmp Rate of change of the internal bubble pressure
+        !!  @param fdMvdt_tmp Rate of change of the mass of vapor in the bubble
         !!  @param advance_EL Rate of change of the mass of vapor in the bubble
-    elemental subroutine s_advance_EL(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, &
-                                      fmass_g, fbeta_c, fbeta_t, fdPbdt_tmp, advance_EL)
+    subroutine s_advance_EL(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, &
+                            fmass_g, fbeta_c, fbeta_t, fdPbdt_tmp, advance_EL)
         $:GPU_ROUTINE(parallelism='[seq]')
-        real(wp), intent(IN) :: fR_tmp, fV_tmp, fPb_tmp, fMv_tmp
-        real(wp), intent(IN) :: fmass_g, fbeta_c, fbeta_t
-        integer, intent(IN) :: bub_id
-        real(wp), intent(INOUT) :: fdPbdt_tmp
+        real(wp), intent(in) :: fR_tmp, fV_tmp, fPb_tmp, fMv_tmp
+        real(wp), intent(in) :: fmass_g, fbeta_c, fbeta_t
+        integer, intent(in) :: bub_id
+        real(wp), intent(inout) :: fdPbdt_tmp
         real(wp), intent(out) :: advance_EL
         real(wp) :: fVapFlux, myR_m, mygamma_m
 

src/simulation/m_bubbles_EE.fpp

@@ -183,7 +183,7 @@ contains
 
         integer :: i, j, k, l, q, ii !< Loop variables
 
-        integer :: adap_dt_stop_max, adap_dt_stop !< Fail-safe exit if max iteration count reached
+        integer :: adap_dt_stop_sum, adap_dt_stop !< Fail-safe exit if max iteration count reached
         integer :: dmBub_id !< Dummy variables for unified subgrid bubble subroutines
         real(wp) :: dmMass_v, dmMass_n, dmBeta_c, dmBeta_t, dmCson
 
@@ -205,10 +205,9 @@ contains
         end do
         $:END_GPU_PARALLEL_LOOP()
 
-        adap_dt_stop_max = 0
+        adap_dt_stop_sum = 0
         $:GPU_PARALLEL_LOOP(private='[j,k,l,Rtmp, Vtmp, myalpha_rho, myalpha, myR, myV, alf, myP, myRho, R2Vav, R3, nbub, pb_local, mv_local, vflux, pbdot, rddot, n_tait, B_tait, my_divu]', collapse=3, &
-            & reduction='[[adap_dt_stop_max]]', reductionOp='[MAX]', &
-            & copy='[adap_dt_stop_max]')
+            & copy='[adap_dt_stop_sum]')
         do l = 0, p
             do k = 0, n
                 do j = 0, m
@@ -298,21 +297,20 @@ contains
                                 pb_local = 0._wp; mv_local = 0._wp; vflux = 0._wp; pbdot = 0._wp
                             end if
 
+                            adap_dt_stop = 0
+
                             ! Adaptive time stepping
                             if (adap_dt) then
-                                adap_dt_stop = 0
 
-                                call s_advance_step(myRho, myP, myR, myV, R0(q), &
-                                                    pb_local, pbdot, alf, n_tait, B_tait, &
-                                                    bub_adv_src(j, k, l), divu_in%sf(j, k, l), &
-                                                    dmBub_id, dmMass_v, dmMass_n, dmBeta_c, &
-                                                    dmBeta_t, dmCson, adap_dt_stop)
+                                adap_dt_stop = f_advance_step(myRho, myP, myR, myV, R0(q), &
+                                                              pb_local, pbdot, alf, n_tait, B_tait, &
+                                                              bub_adv_src(j, k, l), divu_in%sf(j, k, l), &
+                                                              dmBub_id, dmMass_v, dmMass_n, dmBeta_c, &
+                                                              dmBeta_t, dmCson)
 
                                 q_cons_vf(rs(q))%sf(j, k, l) = nbub*myR
                                 q_cons_vf(vs(q))%sf(j, k, l) = nbub*myV
 
-                                adap_dt_stop_max = max(adap_dt_stop_max, adap_dt_stop)
-
                             else
                                 rddot = f_rddot(myRho, myP, myR, myV, R0(q), &
                                                 pb_local, pbdot, alf, n_tait, B_tait, &
@@ -321,14 +319,17 @@ contains
                                 bub_v_src(j, k, l, q) = nbub*rddot
                                 bub_r_src(j, k, l, q) = q_cons_vf(vs(q))%sf(j, k, l)
                             end if
+
+                            $:GPU_ATOMIC(atomic='update')
+                            adap_dt_stop_sum = adap_dt_stop_sum + adap_dt_stop
                         end if
                     end do
                 end do
             end do
         end do
         $:END_GPU_PARALLEL_LOOP()
 
-        if (adap_dt .and. adap_dt_stop_max > 0) call s_mpi_abort("Adaptive time stepping failed to converge.")
+        if (adap_dt .and. adap_dt_stop_sum > 0) call s_mpi_abort("Adaptive time stepping failed to converge.")
 
         if (.not. adap_dt) then
             $:GPU_PARALLEL_LOOP(private='[i,k,l,q]', collapse=3)