diff --git a/Makefile b/Makefile index c2b6800..684c2f9 100755 --- a/Makefile +++ b/Makefile @@ -14,7 +14,7 @@ PROGRAMS = $(MAIN).exe all: $(PROGRAMS) $(PROGRAMS) : $(OBJECTS) - $(FC) $(LDFLAGS) -o $(PROGRAMS) $(OBJECTS) $(LIBFLAGS1) -llapack -lblas $(LIBFLAGS2) -lfftw3 -lm + $(FC) -fopenmp $(LDFLAGS) -o $(PROGRAMS) $(OBJECTS) $(LIBFLAGS1) -llapack -lblas $(LIBFLAGS2) -lfftw3 -lm fftw.o : fftw.f90 $(FC) $(FFLAGS) fftw.f90 @@ -50,7 +50,7 @@ statistics.o : statistics.f90 $(FC) $(FFLAGS) statistics.f90 time_integrators.o : time_integrators.f90 - $(FC) $(FFLAGS) time_integrators.f90 + $(FC) -fopenmp $(FFLAGS) time_integrators.f90 jacobians.o : jacobians.f90 $(FC) $(FFLAGS) jacobians.f90 diff --git a/input.data b/input.data index 2a87028..894b925 100755 --- a/input.data +++ b/input.data @@ -1,7 +1,7 @@ 7.0, 1.5585, 0.1 5.0e+03, 1, 1.1 -20.0, 0.1 +5.0, 0.1 -1.0, 1.0 64, 51, 1 0, 0, 0 -0, 0, 0 +1, 0, 0 diff --git a/time_integrators.f90 b/time_integrators.f90 old mode 100755 new mode 100644 index ac1b6df..8149103 --- a/time_integrators.f90 +++ b/time_integrators.f90 @@ -6,6 +6,7 @@ module time_integrators use allocate_vars use bc_setup use statistics +use omp_lib integer :: it, jt, kkt integer :: info @@ -31,19 +32,23 @@ subroutine imex_rk(vtk_print, save_nusselt) integer :: nprint logical :: wvtk real(dp) :: nusselt_num - +real(dp) :: start, finish +real(dp) :: start_overall, finish_overall +integer :: nthreads, myid +integer, EXTERNAL :: OMP_GET_THREAD_NUM, OMP_GET_NUM_THREADS +real(dp), EXTERNAL :: OMP_GET_WTIME !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (present(vtk_print)) then - wvtk = .true. - nprint = vtk_print + wvtk = .true. + nprint = vtk_print else - wvtk = .false. - nprint = 100 + wvtk = .false. + nprint = 100 end if if (wvtk) then - call write_to_vtk(0, .false.) ! false = Fourier space + call write_to_vtk(0, .false.) ! false = Fourier space end if dt = dt_init @@ -66,31 +71,36 @@ subroutine imex_rk(vtk_print, save_nusselt) ! Time integration do ! while (time < t_final) + start_overall = OMP_GET_WTIME() - dt_final = t_final - time + dt_final = t_final - time - if (dt_final <= dt) then + if (dt_final <= dt) then time = t_final - else + else time = time + dt - end if - - write(*,*) "time = ", time, "dt = ", dt - - nti = nti + 1 - - !::::::::::: - ! STAGE 1 :: - !::::::::::: - phii = phi - Ti = T - uxi = ux - uyi = uy - call calc_explicit(1) - do it = 1,Nx ! kx loop + end if + + write(*,*) "time = ", time, "dt = ", dt + + nti = nti + 1 + + !::::::::::: + ! STAGE 1 :: + !::::::::::: + phii = phi + Ti = T + uxi = ux + uyi = uy + start = OMP_GET_WTIME() + call calc_explicit(1) + finish = OMP_GET_WTIME() + write(*,*) " - calc_explicit(1) timing: ", finish-start, "(s)" + start = OMP_GET_WTIME() + do it = 1,Nx ! kx loop ! Compute phi1 and T1 call calc_vari(tmp_phi, tmp_T, acoeffs(1,1), 1) - ! Compute v1 from phi1 + ! Compute v1 from phi1 call calc_vi(tmp_uy, tmp_phi) ! BOUNDAY CONDITIONS! call update_bcs(tmp_phi1,tmp_uy1, tmp_phi,tmp_uy) @@ -102,25 +112,31 @@ subroutine imex_rk(vtk_print, save_nusselt) K1_T(:,it) = tmp_K_T ! Compute u1 from v1 if (kx(it) /= 0.0_dp) then - !uxi(:,it) = -CI*d1y(tmp_uy)/kx(it) - uxi(:,it) = CI*d1y(tmp_uy)/kx(it) + !uxi(:,it) = -CI*d1y(tmp_uy)/kx(it) + uxi(:,it) = CI*d1y(tmp_uy)/kx(it) else if (kx(it) == 0.0_dp) then - uxi(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! + uxi(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! end if phii(:,it) = tmp_phi Ti (:,it) = tmp_T uyi (:,it) = tmp_uy - end do - ! Compute K2hat - call calc_explicit(2) - - !::::::::::: - ! STAGE 2 :: - !::::::::::: - do it = 1,Nx ! kx loop + end do + finish = OMP_GET_WTIME() + write(*,*) " - stage 1 mid timing: ", finish-start, "(s)" + ! Compute K2hat + start = OMP_GET_WTIME() + call calc_explicit(2) + finish = OMP_GET_WTIME() + write(*,*) " - calc_explicit(2) timing: ", finish-start, "(s)" + + !::::::::::: + ! STAGE 2 :: + !::::::::::: + start = OMP_GET_WTIME() + do it = 1,Nx ! kx loop ! Compute phi2 and T2 call calc_vari(tmp_phi, tmp_T, acoeffs(2,2), 2) - ! Compute v1 from phi1 + ! Compute v1 from phi1 call calc_vi(tmp_uy, tmp_phi) ! BOUNDAY CONDITIONS! call update_bcs(tmp_phi1,tmp_uy1, tmp_phi,tmp_uy) @@ -132,22 +148,28 @@ subroutine imex_rk(vtk_print, save_nusselt) K2_T(:,it) = tmp_K_T ! Compute u1 from v1 if (kx(it) /= 0.0_dp) then - !uxi(:,it) = -CI*d1y(tmp_uy)/kx(it) - uxi(:,it) = CI*d1y(tmp_uy)/kx(it) + !uxi(:,it) = -CI*d1y(tmp_uy)/kx(it) + uxi(:,it) = CI*d1y(tmp_uy)/kx(it) else if (kx(it) == 0.0_dp) then - uxi(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! + uxi(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! end if phii(:,it) = tmp_phi Ti (:,it) = tmp_T uyi (:,it) = tmp_uy - end do - ! Compute K3hat - call calc_explicit(3) - - !::::::::::: - ! STAGE 3 :: - !::::::::::: - do it = 1,Nx ! kx loop + end do + finish = OMP_GET_WTIME() + write(*,*) " - stage 2 mid timing: ", finish-start, "(s)" + ! Compute K3hat + start = OMP_GET_WTIME() + call calc_explicit(3) + finish = OMP_GET_WTIME() + write(*,*) " - calc_explicit(3) timing: ", finish-start, "(s)" + + !::::::::::: + ! STAGE 3 :: + !::::::::::: + start = OMP_GET_WTIME() + do it = 1,Nx ! kx loop ! Compute phi3 and T3 call calc_vari(tmp_phi, tmp_T, acoeffs(3,3), 3) ! Compute v1 from phi1 @@ -162,62 +184,75 @@ subroutine imex_rk(vtk_print, save_nusselt) K3_T(:,it) = tmp_K_T ! Compute u1 from v1 if (kx(it) /= 0.0_dp) then - !uxi(:,it) = -CI*d1y(tmp_uy)/kx(it) - uxi(:,it) = CI*d1y(tmp_uy)/kx(it) + !uxi(:,it) = -CI*d1y(tmp_uy)/kx(it) + uxi(:,it) = CI*d1y(tmp_uy)/kx(it) else if (kx(it) == 0.0_dp) then - uxi(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! + uxi(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! end if phii(:,it) = tmp_phi Ti (:,it) = tmp_T uyi (:,it) = tmp_uy - end do - ! Compute K4hat - call calc_explicit(4) - - ! UPDATE SOLUTIONS - - ! Get phi - phi(2:Ny-1,:) = phi(2:Ny-1,:) + dt*(b(1)*(K1_phi(2:Ny-1,:) + K2hat_phi(2:Ny-1,:)) + & + end do + finish = OMP_GET_WTIME() + write(*,*) " - stage 3 mid timing: ", finish-start, "(s)" + ! Compute K4hat + start = OMP_GET_WTIME() + call calc_explicit(4) + finish = OMP_GET_WTIME() + write(*,*) " - calc_explicit(4) timing: ", finish-start, "(s)" + + ! UPDATE SOLUTIONS + + start = OMP_GET_WTIME() + ! Get phi + phi(2:Ny-1,:) = phi(2:Ny-1,:) + dt*(b(1)*(K1_phi(2:Ny-1,:) + K2hat_phi(2:Ny-1,:)) + & & b(2)*(K2_phi(2:Ny-1,:) + K3hat_phi(2:Ny-1,:)) + & & b(3)*(K3_phi(2:Ny-1,:) + K4hat_phi(2:Ny-1,:))) - ! Get temperature - T(2:Ny-1,:) = T(2:Ny-1,:) + dt*(b(1)*(K1_T(2:Ny-1,:) + K2hat_T(2:Ny-1,:)) + & + ! Get temperature + T(2:Ny-1,:) = T(2:Ny-1,:) + dt*(b(1)*(K1_T(2:Ny-1,:) + K2hat_T(2:Ny-1,:)) + & & b(2)*(K2_T(2:Ny-1,:) + K3hat_T(2:Ny-1,:)) + & & b(3)*(K3_T(2:Ny-1,:) + K4hat_T(2:Ny-1,:))) - ! Get ux and uy - do it = 1,Nx + ! Get ux and uy + do it = 1,Nx ! Solve for v call calc_vi(tmp_uy, phi(:,it)) uy(:,it) = tmp_uy ! Solve for u if (kx(it) /= 0.0_dp) then - !ux(:,it) = -CI*d1y(tmp_uy)/kx(it) - ux(:,it) = CI*d1y(tmp_uy)/kx(it) + !ux(:,it) = -CI*d1y(tmp_uy)/kx(it) + ux(:,it) = CI*d1y(tmp_uy)/kx(it) else if (kx(it) == 0.0_dp) then - ux(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! + ux(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! end if - end do + end do + finish = OMP_GET_WTIME() + write(*,*) " - update sols timing: ", finish-start, "(s)" - if (time == t_final) then + if (time == t_final) then exit - end if - - !call update_dt + end if - if (wvtk) then + !call update_dt + if (wvtk) then if (mod(nti,vtk_print) == 0) then - call write_to_vtk(nti, .false.) ! false = Fourier space + call write_to_vtk(nti, .false.) ! false = Fourier space end if - end if + end if - ! Calculate nusselt number. - if (save_nusselt) then + start = OMP_GET_WTIME() + ! Calculate nusselt number. + if (save_nusselt) then call nusselt(nusselt_num, .true.) ! true = Fourier space write(8000, fmt=1000) nusselt_num flush(8000) - end if + end if + finish = OMP_GET_WTIME() + write(*,*) " - nusselt calc timing: ", finish-start, "(s)" + + finish_overall = OMP_GET_WTIME() + write(*,*) "overall timing: ", finish_overall-start_overall, "(s)" end do ! time loop @@ -260,18 +295,18 @@ subroutine calc_vari(phiout,Tout, aii, stage) ! LHS Matrix (tridiagonal, not necessarily symmetric) do jt = 2,Ny-1 - ddT (jt-1) = 1.0_dp - kappa0*dt*aii*(-kx(it)**2.0_dp + g2(jt)) - dphi(jt-1) = 1.0_dp - nu0 *dt*aii*(-kx(it)**2.0_dp + g2(jt)) + ddT (jt-1) = 1.0_dp - kappa0*dt*aii*(-kx(it)**2.0_dp + g2(jt)) + dphi(jt-1) = 1.0_dp - nu0 *dt*aii*(-kx(it)**2.0_dp + g2(jt)) end do do jt = 2,Ny-2 - duT (jt-1) = -kappa0*dt*g3(jt)*aii - duphi(jt-1) = -nu0 *dt*g3(jt)*aii + duT (jt-1) = -kappa0*dt*g3(jt)*aii + duphi(jt-1) = -nu0 *dt*g3(jt)*aii end do do jt = 3,Ny-1 - dlT (jt-2) = -kappa0*dt*g1(jt)*aii - dlphi(jt-2) = -nu0 *dt*g1(jt)*aii + dlT (jt-2) = -kappa0*dt*g1(jt)*aii + dlphi(jt-2) = -nu0 *dt*g1(jt)*aii end do select case (stage) @@ -289,7 +324,7 @@ subroutine calc_vari(phiout,Tout, aii, stage) call dgtsv(Ny-2, 2, dlT, ddT, duT, T_rhs, Ny-2, info) Tout(2:Ny-1) = cmplx(T_rhs(:,1), T_rhs(:,2), kind=C_DOUBLE_COMPLEX) - ! Set temperature boundary conditions + ! Set temperature boundary conditions Tout(1) = T(1,it) Tout(Ny) = T(Ny,it) @@ -314,7 +349,7 @@ subroutine calc_vari(phiout,Tout, aii, stage) call dgtsv(Ny-2, 2, dlT, ddT, duT, T_rhs, Ny-2, info) Tout(2:Ny-1) = cmplx(T_rhs(:,1), T_rhs(:,2), kind=C_DOUBLE_COMPLEX) - ! Set temperature boundary conditions + ! Set temperature boundary conditions Tout(1) = T(1,it) Tout(Ny) = T(Ny,it) @@ -323,27 +358,27 @@ subroutine calc_vari(phiout,Tout, aii, stage) case(3) Fphi = phi(2:Ny-1,it) + dt*(acoeffs(3,1)*K1_phi(2:Ny-1,it) + & - &acoeffs(3,2)*K2_phi(2:Ny-1,it) + & - &ahatcoeffs(4,1)*K1hat_phi(2:Ny-1,it) + & - &ahatcoeffs(4,2)*K2hat_phi(2:Ny-1,it) + & - &ahatcoeffs(4,3)*K3hat_phi(2:Ny-1,it)) + &acoeffs(3,2)*K2_phi(2:Ny-1,it) + & + &ahatcoeffs(4,1)*K1hat_phi(2:Ny-1,it) + & + &ahatcoeffs(4,2)*K2hat_phi(2:Ny-1,it) + & + &ahatcoeffs(4,3)*K3hat_phi(2:Ny-1,it)) FT = T (2:Ny-1,it) + dt*(acoeffs(3,1)*K1_T (2:Ny-1,it) + & - &acoeffs(3,2)*K2_T (2:Ny-1,it) + & - &ahatcoeffs(4,1)*K1hat_T (2:Ny-1,it) + & - &ahatcoeffs(4,2)*K2hat_T (2:Ny-1,it) + & - &ahatcoeffs(4,3)*K3hat_T (2:Ny-1,it)) + &acoeffs(3,2)*K2_T (2:Ny-1,it) + & + &ahatcoeffs(4,1)*K1hat_T (2:Ny-1,it) + & + &ahatcoeffs(4,2)*K2hat_T (2:Ny-1,it) + & + &ahatcoeffs(4,3)*K3hat_T (2:Ny-1,it)) FT(1) = FT(1) + kappa0*dt*aii*g1(2)*T(1,it) ! b/c Ti(y_1) = T(y_1) FT(Ny-2) = FT(Ny-2) + kappa0*dt*aii*g3(Ny-1)*T(Ny,it) ! b/c Ti(Ny) = T(Ny) phi_rhs(:,1) = real(Fphi) phi_rhs(:,2) = aimag(Fphi) - + T_rhs (:,1) = real(FT) T_rhs (:,2) = aimag(FT) call dgtsv(Ny-2, 2, dlT, ddT, duT, T_rhs, Ny-2, info) Tout(2:Ny-1) = cmplx(T_rhs(:,1), T_rhs(:,2), kind=C_DOUBLE_COMPLEX) - ! Set temperature boundary conditions + ! Set temperature boundary conditions Tout(1) = T(1,it) Tout(Ny) = T(Ny,it) @@ -375,111 +410,132 @@ subroutine calc_explicit(stage) integer :: i, j integer, intent(in) :: stage +real(dp) :: start, finish +start = OMP_GET_WTIME() select case(stage) - case (1) + case (1) do i = 1,Nx - K1hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) + K1hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) end do - case (2) + case (2) do i = 1,Nx - K2hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) + K2hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) end do - case (3) + case (3) do i = 1,Nx - K3hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) + K3hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) end do - case (4) + case (4) do i = 1,Nx - K4hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) + K4hat_phi(:,i) = -kx(i)**2.0_dp*Ti(:,i) end do end select +finish = OMP_GET_WTIME() +write(*,*) " - - l1 timing: ", finish-start, "(s)" +start = OMP_GET_WTIME() do i=1,Nx - ! Compute dx(T) in Fourier space - nlT (:,i) = kx(i)*Ti(:,i) - ! Compute D2(ux) - nlphi(:,i) = -kx(i)**2.0_dp*uxi(:,i) + d2y(uxi(:,i)) + ! Compute dx(T) in Fourier space + nlT (:,i) = kx(i)*Ti(:,i) + ! Compute D2(ux) + nlphi(:,i) = -kx(i)**2.0_dp*uxi(:,i) + d2y(uxi(:,i)) end do +finish = OMP_GET_WTIME() +write(*,*) " - - l2 timing: ", finish-start, "(s)" + !nlT = -CI*nlT nlT = CI*nlT +start = OMP_GET_WTIME() do j = 1,Ny - ! Bring everything to physical space - tnlT = nlT(j,:) - tnlphi = nlphi(j,:) - tT = Ti(j,:) - tux = uxi(j,:) - tuy = uyi(j,:) - tphi = phii(j,:) - call fftw_execute_dft(iplannlT, tnlT, tnlT) - call fftw_execute_dft(iplannlphi, tnlphi, tnlphi) - call fftw_execute_dft(iplanT, tT, tT) - call fftw_execute_dft(iplanux, tux, tux) - call fftw_execute_dft(iplanuy, tuy, tuy) - call fftw_execute_dft(iplanphi, tphi, tphi) - nlT(j,:) = tnlT - nlphi(j,:) = tnlphi - Ti(j,:) = tT - uxi(j,:) = tux - uyi(j,:) = tuy - phii(j,:) = tphi + ! Bring everything to physical space + tnlT = nlT(j,:) + tnlphi = nlphi(j,:) + tT = Ti(j,:) + tux = uxi(j,:) + tuy = uyi(j,:) + tphi = phii(j,:) + call fftw_execute_dft(iplannlT, tnlT, tnlT) + call fftw_execute_dft(iplannlphi, tnlphi, tnlphi) + call fftw_execute_dft(iplanT, tT, tT) + call fftw_execute_dft(iplanux, tux, tux) + call fftw_execute_dft(iplanuy, tuy, tuy) + call fftw_execute_dft(iplanphi, tphi, tphi) + nlT(j,:) = tnlT + nlphi(j,:) = tnlphi + Ti(j,:) = tT + uxi(j,:) = tux + uyi(j,:) = tuy + phii(j,:) = tphi end do +finish = OMP_GET_WTIME() +write(*,*) " - - l3 timing: ", finish-start, "(s)" ! Calculate nonlinear term +start = OMP_GET_WTIME() do i = 1,Nx - ! Temperature - tmp_T = Ti(:,i) - nlT(:,i) = uxi(:,i)*nlT(:,i) + uyi(:,i)*d1y(tmp_T) - ! phi - nlphi(:,i) = uxi(:,i)*phii(:,i) - uyi(:,i)*nlphi(:,i) + ! Temperature + tmp_T = Ti(:,i) + nlT(:,i) = uxi(:,i)*nlT(:,i) + uyi(:,i)*d1y(tmp_T) + ! phi + nlphi(:,i) = uxi(:,i)*phii(:,i) - uyi(:,i)*nlphi(:,i) end do +finish = OMP_GET_WTIME() +write(*,*) " - - l4 timing: ", finish-start, "(s)" ! Bring nonlinear terms back to Fourier space +start = OMP_GET_WTIME() do j = 1,Ny - tnlT = nlT(j,:) - tnlphi = nlphi(j,:) - call fftw_execute_dft(plannlT, tnlT, tnlT) - call fftw_execute_dft(plannlphi, tnlphi, tnlphi) - ! Dealias - do i = 1,Nx + tnlT = nlT(j,:) + tnlphi = nlphi(j,:) + call fftw_execute_dft(plannlT, tnlT, tnlT) + call fftw_execute_dft(plannlphi, tnlphi, tnlphi) + ! Dealias + do i = 1,Nx if (abs(kx(i))/alpha >= Nf/2) then - tnlT(i) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) - tnlphi(i) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) + tnlT(i) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) + tnlphi(i) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) end if - end do - nlT(j,:) = tnlT - nlphi(j,:) = tnlphi + end do + nlT(j,:) = tnlT + nlphi(j,:) = tnlphi end do +finish = OMP_GET_WTIME() +write(*,*) " - - l5 timing: ", finish-start, "(s)" + nlT = nlT / real(Nx,kind=dp) nlphi = nlphi / real(Nx,kind=dp) +start = OMP_GET_WTIME() select case (stage) - case (1) + case (1) do i = 1,Nx !K1hat_phi(:,i) = K1hat_phi(:,i) + CI*kx(i)*nlphi(:,i) K1hat_phi(:,i) = K1hat_phi(:,i) - CI*kx(i)*nlphi(:,i) end do K1hat_T = -nlT - case (2) + case (2) do i = 1,Nx !K2hat_phi(:,i) = K2hat_phi(:,i) + CI*kx(i)*nlphi(:,i) K2hat_phi(:,i) = K2hat_phi(:,i) - CI*kx(i)*nlphi(:,i) end do K2hat_T = -nlT - case (3) + case (3) do i = 1,Nx !K3hat_phi(:,i) = K3hat_phi(:,i) + CI*kx(i)*nlphi(:,i) K3hat_phi(:,i) = K3hat_phi(:,i) - CI*kx(i)*nlphi(:,i) end do K3hat_T = -nlT - case (4) + case (4) do i = 1,Nx !K4hat_phi(:,i) = K4hat_phi(:,i) + CI*kx(i)*nlphi(:,i) K4hat_phi(:,i) = K4hat_phi(:,i) - CI*kx(i)*nlphi(:,i) end do K4hat_T = -nlT end select +finish = OMP_GET_WTIME() +write(*,*) " - - l6 timing: ", finish-start, "(s)" end subroutine calc_explicit @@ -508,11 +564,11 @@ subroutine calc_vi(vi, phiin) duvi = 0.0_dp do j = 2,Ny-1 - dvi(j-1) = -kx(it)**2.0_dp + g2(j) + dvi(j-1) = -kx(it)**2.0_dp + g2(j) end do do j = 2,Ny-2 - duvi(j-1) = g3(j) + duvi(j-1) = g3(j) end do do j = 3,Ny-1 @@ -562,12 +618,12 @@ subroutine update_bcs(phiout,vout, phiin,vin) ! Find c1 and c2. if (detC == 0.0_dp) then - c1 = (0.0_dp, 0.0_dp) - c2 = (0.0_dp, 0.0_dp) + c1 = (0.0_dp, 0.0_dp) + c2 = (0.0_dp, 0.0_dp) else - c1t = (C(2,2)*c1 - C(1,2)*c2) / detC - c2 = (C(1,1)*c2 - C(2,1)*c1) / detC - c1 = c1t + c1t = (C(2,2)*c1 - C(1,2)*c2) / detC + c2 = (C(1,1)*c2 - C(2,1)*c1) / detC + c1 = c1t end if ! Update uy and Phi. @@ -586,38 +642,38 @@ subroutine update_dt uyi = uy do jj = 1,Ny - ! Bring everything to physical space - tux = uxi(jj,:) - tuy = uyi(jj,:) - call fftw_execute_dft(iplanux, tux, tux) - call fftw_execute_dft(iplanuy, tuy, tuy) - uxi(jj,:) = tux - uyi(jj,:) = tuy + ! Bring everything to physical space + tux = uxi(jj,:) + tuy = uyi(jj,:) + call fftw_execute_dft(iplanux, tux, tux) + call fftw_execute_dft(iplanuy, tuy, tuy) + uxi(jj,:) = tux + uyi(jj,:) = tuy end do dt_old = dt do jj = 1,Ny - do ii = 1,Nx + do ii = 1,Nx tmp = real(uxi(jj,ii)) / dxmin + real(uyi(jj,ii)) / dymin if (tmp > tmpmax) then - tmpmax = tmp + tmpmax = tmp end if - end do + end do end do dt = cfl / tmpmax if (dt > dt_ramp * dt_old) then - dt = dt_ramp * dt_old + dt = dt_ramp * dt_old else if (dt < dt_old / dt_ramp) then - dt = dt_old / dt_ramp + dt = dt_old / dt_ramp end if if (dt > dtmax) then - dt = dtmax + dt = dtmax else if (dt < dtmin) then - dt = dtmin + dt = dtmin end if call init_bc(acoeffs(1,1))