diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp index 857225007..21f7aaf58 100644 --- a/src/simulation/m_bubbles.fpp +++ b/src/simulation/m_bubbles.fpp @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 5a1927f0c..d897693e2 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -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 diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index ed5628d9c..aace6fdb3 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -590,11 +590,10 @@ contains !> 3rd order TVD RK time-stepping algorithm !! @param t_step Current time-step - subroutine s_3rd_order_tvd_rk(t_step, time_avg, dt_in) + subroutine s_3rd_order_tvd_rk(t_step, time_avg) ! -------------------------------- - integer, intent(in) :: t_step - real(kind(0d0)), intent(inout) :: time_avg - real(kind(0d0)), intent(in) :: dt_in + integer, intent(IN) :: t_step + real(kind(0d0)), intent(INOUT) :: time_avg integer :: i, j, k, l, q !< Generic loop iterator real(kind(0d0)) :: ts_error, denom, error_fraction, time_step_factor !< Generic loop iterator @@ -635,7 +634,7 @@ contains do j = 0, m q_cons_ts(2)%vf(i)%sf(j, k, l) = & q_cons_ts(1)%vf(i)%sf(j, k, l) & - + dt_in*rhs_vf(i)%sf(j, k, l) + + dt*rhs_vf(i)%sf(j, k, l) end do end do end do @@ -651,7 +650,7 @@ contains do q = 1, nnode pb_ts(2)%sf(j, k, l, q, i) = & pb_ts(1)%sf(j, k, l, q, i) & - + dt_in*rhs_pb(j, k, l, q, i) + + dt*rhs_pb(j, k, l, q, i) end do end do end do @@ -668,7 +667,7 @@ contains do q = 1, nnode mv_ts(2)%sf(j, k, l, q, i) = & mv_ts(1)%sf(j, k, l, q, i) & - + dt_in*rhs_mv(j, k, l, q, i) + + dt*rhs_mv(j, k, l, q, i) end do end do end do @@ -710,7 +709,7 @@ contains q_cons_ts(2)%vf(i)%sf(j, k, l) = & (3d0*q_cons_ts(1)%vf(i)%sf(j, k, l) & + q_cons_ts(2)%vf(i)%sf(j, k, l) & - + dt_in*rhs_vf(i)%sf(j, k, l))/4d0 + + dt*rhs_vf(i)%sf(j, k, l))/4d0 end do end do end do @@ -726,7 +725,7 @@ contains pb_ts(2)%sf(j, k, l, q, i) = & (3d0*pb_ts(1)%sf(j, k, l, q, i) & + pb_ts(2)%sf(j, k, l, q, i) & - + dt_in*rhs_pb(j, k, l, q, i))/4d0 + + dt*rhs_pb(j, k, l, q, i))/4d0 end do end do end do @@ -744,7 +743,7 @@ contains mv_ts(2)%sf(j, k, l, q, i) = & (3d0*mv_ts(1)%sf(j, k, l, q, i) & + mv_ts(2)%sf(j, k, l, q, i) & - + dt_in*rhs_mv(j, k, l, q, i))/4d0 + + dt*rhs_mv(j, k, l, q, i))/4d0 end do end do end do @@ -785,7 +784,7 @@ contains q_cons_ts(1)%vf(i)%sf(j, k, l) = & (q_cons_ts(1)%vf(i)%sf(j, k, l) & + 2d0*q_cons_ts(2)%vf(i)%sf(j, k, l) & - + 2d0*dt_in*rhs_vf(i)%sf(j, k, l))/3d0 + + 2d0*dt*rhs_vf(i)%sf(j, k, l))/3d0 end do end do end do @@ -801,7 +800,7 @@ contains pb_ts(1)%sf(j, k, l, q, i) = & (pb_ts(1)%sf(j, k, l, q, i) & + 2d0*pb_ts(2)%sf(j, k, l, q, i) & - + 2d0*dt_in*rhs_pb(j, k, l, q, i))/3d0 + + 2d0*dt*rhs_pb(j, k, l, q, i))/3d0 end do end do end do @@ -819,7 +818,7 @@ contains mv_ts(1)%sf(j, k, l, q, i) = & (mv_ts(1)%sf(j, k, l, q, i) & + 2d0*mv_ts(2)%sf(j, k, l, q, i) & - + 2d0*dt_in*rhs_mv(j, k, l, q, i))/3d0 + + 2d0*dt*rhs_mv(j, k, l, q, i))/3d0 end do end do end do @@ -874,13 +873,13 @@ contains call nvtxStartRange("Time_Step") ! Stage 1 of 3 ===================================================== - call s_3rd_order_tvd_rk(t_step, time_avg, 0.5d0*dt) + call s_adaptive_dt_bubble(t_step) ! Stage 2 of 3 ===================================================== - call s_adaptive_dt_bubble(t_step) + call s_3rd_order_tvd_rk(t_step, time_avg) ! Stage 3 of 3 ===================================================== - call s_3rd_order_tvd_rk(t_step, time_avg, 0.5d0*dt) + call s_adaptive_dt_bubble(t_step) call nvtxEndRange @@ -913,7 +912,9 @@ contains call s_compute_bubble_source(q_cons_ts(1)%vf, q_prim_vf, t_step, rhs_vf) - end subroutine s_adaptive_dt_bubble + call s_comp_alpha_from_n(q_cons_ts(1)%vf) + + end subroutine s_adaptive_dt_bubble ! ------------------------------ !> This subroutine applies the body forces source term at each !! Runge-Kutta stage