Skip to content

Commit

Permalink
Update on adaptive time stepping for sub-grid bubbles (#408)
Browse files Browse the repository at this point in the history 
Co-authored-by: Hyeoksu Lee <hlee3@dt-login04.delta.ncsa.illinois.edu>
  • Loading branch information
@hyeoksu-lee
hyeoksu-lee and Hyeoksu Lee authored Jul 12, 2024
1 parent a9c0e32 commit 78fd0c5
Show file tree
Hide file tree
Showing 3 changed files with 200 additions and 115 deletions.
278 changes: 181 additions & 97 deletions src/simulation/m_bubbles.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,19 +184,13 @@ contains
!! @param q_prim_vf Primitive variables
!! @param q_cons_vf Conservative variables
subroutine s_compute_bubble_source(q_cons_vf, q_prim_vf, t_step, rhs_vf)

type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
integer, intent(in) :: t_step
type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf

real(kind(0d0)) :: tmp1, tmp2, tmp3, tmp4, &
c_gas, c_liquid, &
Cpbw, Cpinf, Cpinf_dot, &
myH, myHdot, rddot, alf_gas

real(kind(0d0)) :: pb, mv, vflux, pldot, pbdot

real(kind(0d0)) :: rddot
real(kind(0d0)) :: pb, mv, vflux, pbdot
real(kind(0d0)) :: n_tait, B_tait

real(kind(0d0)), dimension(nb) :: Rtmp, Vtmp
Expand All @@ -211,11 +205,10 @@ contains
integer :: i, j, k, l, q, ii !< Loop variables
integer :: ndirs !< Number of coordinate directions

real(kind(0d0)) :: err_R, err_V, err !< Error estimates for adaptive time stepping
real(kind(0d0)) :: err1, err2, err3, err4, err5 !< Error estimates for adaptive time stepping
real(kind(0d0)) :: t_new !< Updated time step size
real(kind(0d0)) :: h, h0, h1, h_min !< Time step size
real(kind(0d0)) :: d0, d1, d2 !< norms
real(kind(0d0)), dimension(4) :: myR_tmp, myV_tmp, myA_tmp !< Bubble radius, radial velocity, and radial acceleration for the inner loop
real(kind(0d0)) :: h !< Time step size
real(kind(0d0)), dimension(4) :: myR_tmp1, myV_tmp1, myR_tmp2, myV_tmp2 !< Bubble radius, radial velocity, and radial acceleration for the inner loop

!$acc parallel loop collapse(3) gang vector default(present)
do l = 0, p
Expand All @@ -234,7 +227,7 @@ contains
end do
end do

!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp, myalpha_rho, myalpha, myR_tmp, myV_tmp, myA_tmp)
!$acc parallel loop collapse(3) gang vector default(present) private(Rtmp, Vtmp, myalpha_rho, myalpha, myR_tmp1, myV_tmp1, myR_tmp2, myV_tmp2)
do l = 0, p
do k = 0, n
do j = 0, m
Expand Down Expand Up @@ -327,114 +320,73 @@ contains
! Adaptive time stepping
if (adap_dt) then
! Determine the starting time step
! Evaluate f(x0,y0)
myR_tmp(1) = myR
myV_tmp(1) = myV
myA_tmp(1) = f_rddot(myRho, myP, myR_tmp(1), myV_tmp(1), R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l))

! Compute d0 = ||y0|| and d1 = ||f(x0,y0)||
d0 = dsqrt((myR_tmp(1)**2d0 + myV_tmp(1)**2d0)/2d0)
d1 = dsqrt((myV_tmp(1)**2d0 + myA_tmp(1)**2d0)/2d0)
if (d0 < 1d-5 .or. d1 < 1d-5) then
h0 = 1d-6
else
h0 = 1d-2*(d0/d1)
end if

! Evaluate f(x0+h0,y0+h0*f(x0,y0))
myR_tmp(2) = myR_tmp(1) + h0*myV_tmp(1)
myV_tmp(2) = myV_tmp(1) + h0*myA_tmp(1)
myA_tmp(2) = f_rddot(myRho, myP, myR_tmp(2), myV_tmp(2), R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l))

! Compute d2 = ||f(x0+h0,y0+h0*f(x0,y0))-f(x0,y0)||/h0
d2 = dsqrt(((myV_tmp(2) - myV_tmp(1))**2d0 + (myA_tmp(2) - myA_tmp(1))**2d0)/2d0)/h0

! Set h1 = (0.01/max(d1,d2))^{1/(p+1)}
! if max(d1,d2) < 1e-15, h1 = max(1e-6, h0*1e-3)
if (max(d1, d2) < 1d-15) then
h1 = max(1d-6, h0*1d-3)
else
h1 = (1d-2/max(d1, d2))**(1d0/3d0)
end if

! Set h = min(100*h0,h1)
h = min(100d0*h0, h1)
call s_initialize_adap_dt(myRho, myP, myR, myV, R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l), h)

! Advancing one step
t_new = 0d0
h_min = h
do while (.true.)
if (t_new + h > dt) then
h = dt - t_new
if (t_new + h > 0.5d0*dt) then
h = 0.5d0*dt - t_new
end if

! Advancing one sub-step
do while (.true.)
! Advance one sub-step and evaluate the error
! Stage 0
myR_tmp(1) = myR
myV_tmp(1) = myV
myA_tmp(1) = f_rddot(myRho, myP, myR_tmp(1), myV_tmp(1), R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l))

! Stage 1
myR_tmp(2) = myR_tmp(1) + h*myV_tmp(1)
myV_tmp(2) = myV_tmp(1) + h*myA_tmp(1)
myA_tmp(2) = f_rddot(myRho, myP, myR_tmp(2), myV_tmp(2), R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l))

! Stage 2
myR_tmp(3) = myR_tmp(1) + (h/4d0)*(myV_tmp(1) + myV_tmp(2))
myV_tmp(3) = myV_tmp(1) + (h/4d0)*(myA_tmp(1) + myA_tmp(2))
myA_tmp(3) = f_rddot(myRho, myP, myR_tmp(3), myV_tmp(3), R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l))

! Stage 3
myR_tmp(4) = myR + (h/6d0)*(myV_tmp(1) + myV_tmp(2) + 4*myV_tmp(3))
myV_tmp(4) = myV + (h/6d0)*(myA_tmp(1) + myA_tmp(2) + 4*myA_tmp(3))
myA_tmp(4) = f_rddot(myRho, myP, myR_tmp(4), myV_tmp(4), R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l))

! Estimate error
err_R = (-5d0*h/24d0)*(myV_tmp(2) + myV_tmp(3) - 2d0*myV_tmp(4)) &
/max(abs(myR_tmp(1)), abs(myR_tmp(4)))
err_V = (-5d0*h/24d0)*(myA_tmp(2) + myA_tmp(3) - 2d0*myA_tmp(4)) &
/max(abs(myV_tmp(1)), abs(myV_tmp(4)))
err = dsqrt((err_R**2d0 + err_V**2d0)/2d0)/1d-2 ! Rtol = 1e-2
! Advance one sub-step
call s_advance_substep(myRho, myP, myR, myV, R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l), h, &
myR_tmp1, myV_tmp1, err1)

! Advance one sub-step by advancing two half steps
call s_advance_substep(myRho, myP, myR, myV, R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l), 0.5d0*h, &
myR_tmp2, myV_tmp2, err2)

call s_advance_substep(myRho, myP, myR_tmp2(4), myV_tmp2(4), R0(q), &
pb, pbdot, alf, n_tait, B_tait, &
bub_adv_src(j, k, l), divu%sf(j, k, l), 0.5d0*h, &
myR_tmp2, myV_tmp2, err3)

err4 = abs((myR_tmp1(4) - myR_tmp2(4))/myR_tmp1(4))
err5 = abs((myV_tmp1(4) - myV_tmp2(4))/myV_tmp1(4))
if (abs(myV_tmp1(4)) < 1e-12) err5 = 0d0

! Determine acceptance/rejection and update step size
if ((err <= 1d0) .and. myR_tmp(4) > 0d0) then
! Rule 1: err1, err2, err3 < tol
! Rule 2: myR_tmp1(4) > 0d0
! Rule 3: abs((myR_tmp1(4) - myR_tmp2(4))/myR) < tol
! Rule 4: abs((myV_tmp1(4) - myV_tmp2(4))/myV) < tol
if ((err1 <= 1d-4) .and. (err2 <= 1d-4) .and. (err3 <= 1d-4) &
.and. (err4 < 1d-4) .and. (err5 < 1d-4) &
.and. myR_tmp1(4) > 0d0) then

! Accepted. Finalize the sub-step
myR = myR_tmp(4)
myV = myV_tmp(4)
t_new = t_new + h

! Update R and V
myR = myR_tmp1(4)
myV = myV_tmp1(4)

! Update step size for the next sub-step
h = h*min(2d0, max(0.5d0, (1d0/err)**(1d0/3d0)))
h = h*min(2d0, max(0.5d0, (1d-4/err1)**(1d0/3d0)))

exit
else
! Rejected. Update step size for the next try on sub-step
if (myR_tmp(4) > 0d0) then
h = h*min(2d0, max(0.5d0, (1d0/err)**(1d0/3d0)))
else
if (err2 <= 1d-4) then
h = 0.5d0*h
else
h = 0.25d0*h
end if
end if

if (h < h_min) h_min = h
end if
end do

! Exit the loop if the final time reached dt
if (t_new == dt) exit
if (t_new == 0.5d0*dt) exit

end do

Expand Down Expand Up @@ -486,6 +438,138 @@ contains
end if
end subroutine s_compute_bubble_source

!> Choose the initial time step size for the adaptive time stepping routine
!! (See Heirer, E. Hairer S.P.Nørsett G. Wanner, Solving Ordinary
!! Differential Equations I, Chapter II.4)
!! @param fRho Current density
!! @param fP Current driving pressure
!! @param fR Current bubble radius
!! @param fV Current bubble velocity
!! @param fR0 Equilibrium bubble radius
!! @param fpb Internal bubble pressure
!! @param fpbdot Time-derivative of internal bubble pressure
!! @param alf bubble volume fraction
!! @param fntait Tait EOS parameter
!! @param fBtait Tait EOS parameter
!! @param f_bub_adv_src Source for bubble volume fraction
!! @param f_divu Divergence of velocity
!! @param h Time step size
subroutine s_initialize_adap_dt(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
fntait, fBtait, f_bub_adv_src, f_divu, h)
!$acc routine seq
real(kind(0d0)), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
real(kind(0d0)), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu
real(kind(0d0)), intent(out) :: h

real(kind(0d0)) :: h0, h1, h_min !< Time step size
real(kind(0d0)) :: d0, d1, d2 !< norms
real(kind(0d0)), dimension(2) :: myR_tmp, myV_tmp, myA_tmp !< Bubble radius, radial velocity, and radial acceleration

! Determine the starting time step
! Evaluate f(x0,y0)
myR_tmp(1) = fR
myV_tmp(1) = fV
myA_tmp(1) = f_rddot(fRho, fP, myR_tmp(1), myV_tmp(1), fR0, &
fpb, fpbdot, alf, fntait, fBtait, &
f_bub_adv_src, f_divu)

! Compute d0 = ||y0|| and d1 = ||f(x0,y0)||
d0 = DSQRT((myR_tmp(1)**2d0 + myV_tmp(1)**2d0)/2d0)
d1 = DSQRT((myV_tmp(1)**2d0 + myA_tmp(1)**2d0)/2d0)
if (d0 < 1d-5 .or. d1 < 1d-5) then
h0 = 1d-6
else
h0 = 1d-2*(d0/d1)
end if

! Evaluate f(x0+h0,y0+h0*f(x0,y0))
myR_tmp(2) = myR_tmp(1) + h0*myV_tmp(1)
myV_tmp(2) = myV_tmp(1) + h0*myA_tmp(1)
myA_tmp(2) = f_rddot(fRho, fP, myR_tmp(2), myV_tmp(2), fR0, &
fpb, fpbdot, alf, fntait, fBtait, &
f_bub_adv_src, f_divu)

! Compute d2 = ||f(x0+h0,y0+h0*f(x0,y0))-f(x0,y0)||/h0
d2 = DSQRT(((myV_tmp(2) - myV_tmp(1))**2d0 + (myA_tmp(2) - myA_tmp(1))**2d0)/2d0)/h0

! Set h1 = (0.01/max(d1,d2))^{1/(p+1)}
! if max(d1,d2) < 1e-15, h1 = max(1e-6, h0*1e-3)
if (max(d1, d2) < 1d-15) then
h1 = max(1d-6, h0*1d-3)
else
h1 = (1d-2/max(d1, d2))**(1d0/3d0)
end if

! Set h = min(100*h0,h1)
h = min(100d0*h0, h1)

end subroutine s_initialize_adap_dt

!> Integrate bubble variables over the given time step size, h
!! @param fRho Current density
!! @param fP Current driving pressure
!! @param fR Current bubble radius
!! @param fV Current bubble velocity
!! @param fR0 Equilibrium bubble radius
!! @param fpb Internal bubble pressure
!! @param fpbdot Time-derivative of internal bubble pressure
!! @param alf bubble volume fraction
!! @param fntait Tait EOS parameter
!! @param fBtait Tait EOS parameter
!! @param f_bub_adv_src Source for bubble volume fraction
!! @param f_divu Divergence of velocity
!! @param h Time step size
!! @param myR_tmp Bubble radius at each stage
!! @param myV_tmp Bubble radial velocity at each stage
!! @param err Estimated error
subroutine s_advance_substep(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
fntait, fBtait, f_bub_adv_src, f_divu, h, &
myR_tmp, myV_tmp, err)
!$acc routine seq
real(kind(0d0)), intent(IN) :: fRho, fP, fR, fV, fR0, fpb, fpbdot, alf
real(kind(0d0)), intent(IN) :: fntait, fBtait, f_bub_adv_src, f_divu, h
real(kind(0d0)), dimension(4), intent(OUT) :: myR_tmp, myV_tmp
real(kind(0d0)), dimension(4) :: myA_tmp
real(kind(0d0)), intent(OUT) :: err
real(kind(0d0)) :: err_R, err_V

! Stage 0
myR_tmp(1) = fR
myV_tmp(1) = fV
myA_tmp(1) = f_rddot(fRho, fP, myR_tmp(1), myV_tmp(1), fR0, &
fpb, fpbdot, alf, fntait, fBtait, &
f_bub_adv_src, f_divu)

! Stage 1
myR_tmp(2) = myR_tmp(1) + h*myV_tmp(1)
myV_tmp(2) = myV_tmp(1) + h*myA_tmp(1)
myA_tmp(2) = f_rddot(fRho, fP, myR_tmp(2), myV_tmp(2), fR0, &
fpb, fpbdot, alf, fntait, fBtait, &
f_bub_adv_src, f_divu)

! Stage 2
myR_tmp(3) = myR_tmp(1) + (h/4d0)*(myV_tmp(1) + myV_tmp(2))
myV_tmp(3) = myV_tmp(1) + (h/4d0)*(myA_tmp(1) + myA_tmp(2))
myA_tmp(3) = f_rddot(fRho, fP, myR_tmp(3), myV_tmp(3), fR0, &
fpb, fpbdot, alf, fntait, fBtait, &
f_bub_adv_src, f_divu)

! Stage 3
myR_tmp(4) = myR_tmp(1) + (h/6d0)*(myV_tmp(1) + myV_tmp(2) + 4d0*myV_tmp(3))
myV_tmp(4) = myV_tmp(1) + (h/6d0)*(myA_tmp(1) + myA_tmp(2) + 4d0*myA_tmp(3))
myA_tmp(4) = f_rddot(fRho, fP, myR_tmp(4), myV_tmp(4), fR0, &
fpb, fpbdot, alf, fntait, fBtait, &
f_bub_adv_src, f_divu)

! Estimate error
err_R = (-5d0*h/24d0)*(myV_tmp(2) + myV_tmp(3) - 2d0*myV_tmp(4)) &
/max(abs(myR_tmp(1)), abs(myR_tmp(4)))
err_V = (-5d0*h/24d0)*(myA_tmp(2) + myA_tmp(3) - 2d0*myA_tmp(4)) &
/max(abs(myV_tmp(1)), abs(myV_tmp(4)))
err = DSQRT((err_R**2d0 + err_V**2d0)/2d0)

end subroutine s_advance_substep

!> Function that computes that bubble wall pressure for Gilmore bubbles
!! @param fR0 Equilibrium bubble radius
!! @param fR Current bubble radius
Expand Down
2 changes: 1 addition & 1 deletion src/simulation/m_start_up.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1113,7 +1113,7 @@ contains
elseif (time_stepper == 2) then
call s_2nd_order_tvd_rk(t_step, time_avg)
elseif (time_stepper == 3 .and. (.not. adap_dt)) then
call s_3rd_order_tvd_rk(t_step, time_avg, dt)
call s_3rd_order_tvd_rk(t_step, time_avg)
elseif (time_stepper == 3 .and. adap_dt) then
call s_strang_splitting(t_step, time_avg)
end if
Expand Down
Loading

0 comments on commit 78fd0c5

Please sign in to comment.