From b363234e419557fc418724c99339e892ebf2b7a4 Mon Sep 17 00:00:00 2001 From: student 328724 group 84102 class g_84102 Date: Tue, 13 Apr 2021 16:25:17 -0400 Subject: [PATCH 1/9] First commit from the academic cluster. Adding open_mp timing for calc explicit calls and entire loop. Timing results show that 60% of runtime comes from calc_explicit call --- Makefile | 4 ++-- input.data | 4 ++-- time_integrators.f90 | 26 ++++++++++++++++++++++++-- 3 files changed, 28 insertions(+), 6 deletions(-) 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..e04fbfc 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 +100.0, 0.1 -1.0, 1.0 -64, 51, 1 +640, 510, 1 0, 0, 0 0, 0, 0 diff --git a/time_integrators.f90 b/time_integrators.f90 index ac1b6df..c4e9d40 100755 --- 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,7 +32,11 @@ 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 @@ -66,6 +71,7 @@ subroutine imex_rk(vtk_print, save_nusselt) ! Time integration do ! while (time < t_final) + start_overall = OMP_GET_WTIME() dt_final = t_final - time @@ -86,7 +92,10 @@ subroutine imex_rk(vtk_print, save_nusselt) 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)" do it = 1,Nx ! kx loop ! Compute phi1 and T1 call calc_vari(tmp_phi, tmp_T, acoeffs(1,1), 1) @@ -112,7 +121,10 @@ subroutine imex_rk(vtk_print, save_nusselt) uyi (:,it) = tmp_uy end do ! Compute K2hat + start = OMP_GET_WTIME() call calc_explicit(2) + finish = OMP_GET_WTIME() + write(*,*) " - calc_explicit(2) timing: ", finish-start, "(s)" !::::::::::: ! STAGE 2 :: @@ -142,7 +154,10 @@ subroutine imex_rk(vtk_print, save_nusselt) uyi (:,it) = tmp_uy end do ! Compute K3hat + start = OMP_GET_WTIME() call calc_explicit(3) + finish = OMP_GET_WTIME() + write(*,*) " - calc_explicit(3) timing: ", finish-start, "(s)" !::::::::::: ! STAGE 3 :: @@ -172,7 +187,10 @@ subroutine imex_rk(vtk_print, save_nusselt) uyi (:,it) = tmp_uy end do ! Compute K4hat + start = OMP_GET_WTIME() call calc_explicit(4) + finish = OMP_GET_WTIME() + write(*,*) " - calc_explicit(4) timing: ", finish-start, "(s)" ! UPDATE SOLUTIONS @@ -205,7 +223,8 @@ subroutine imex_rk(vtk_print, save_nusselt) end if !call update_dt - + ! Don't write vtk for now. + wvtk = .false. if (wvtk) then if (mod(nti,vtk_print) == 0) then call write_to_vtk(nti, .false.) ! false = Fourier space @@ -218,6 +237,9 @@ subroutine imex_rk(vtk_print, save_nusselt) write(8000, fmt=1000) nusselt_num flush(8000) end if + + finish_overall = OMP_GET_WTIME() + write(*,*) "overall timing: ", finish_overall-start_overall, "(s)" end do ! time loop From 5118c5cde2b591abc1af62e68c09e838c1e93090 Mon Sep 17 00:00:00 2001 From: student 328724 group 84102 class g_84102 Date: Tue, 13 Apr 2021 16:57:10 -0400 Subject: [PATCH 2/9] Adding more refined timing data to include loops inside the calc_explicit function --- input.data | 2 +- time_integrators.f90 | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/input.data b/input.data index e04fbfc..36d6ad1 100755 --- a/input.data +++ b/input.data @@ -2,6 +2,6 @@ 5.0e+03, 1, 1.1 100.0, 0.1 -1.0, 1.0 -640, 510, 1 +1600, 1275, 1 0, 0, 0 0, 0, 0 diff --git a/time_integrators.f90 b/time_integrators.f90 index c4e9d40..1f966e9 100755 --- a/time_integrators.f90 +++ b/time_integrators.f90 @@ -397,7 +397,9 @@ subroutine calc_explicit(stage) integer :: i, j integer, intent(in) :: stage +real(dp) :: start, finish +start = OMP_GET_WTIME() select case(stage) case (1) do i = 1,Nx @@ -416,16 +418,23 @@ subroutine calc_explicit(stage) 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)) 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,:) @@ -447,8 +456,11 @@ subroutine calc_explicit(stage) 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) @@ -456,8 +468,11 @@ subroutine calc_explicit(stage) ! 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,:) @@ -473,9 +488,13 @@ subroutine calc_explicit(stage) 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) do i = 1,Nx @@ -502,6 +521,8 @@ subroutine calc_explicit(stage) end do K4hat_T = -nlT end select +finish = OMP_GET_WTIME() +write(*,*) " - - l6 timing: ", finish-start, "(s)" end subroutine calc_explicit From 2c1cd482f6a7588827e98ff8716495cec80e2919 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Fri, 16 Apr 2021 21:30:22 +0000 Subject: [PATCH 3/9] Modify input --- input.data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/input.data b/input.data index 36d6ad1..ca77d7c 100755 --- a/input.data +++ b/input.data @@ -1,6 +1,6 @@ 7.0, 1.5585, 0.1 5.0e+03, 1, 1.1 -100.0, 0.1 +5.0, 0.1 -1.0, 1.0 1600, 1275, 1 0, 0, 0 From e87066a51a5a5fc526ccf414123144dbd82a5d3f Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sun, 18 Apr 2021 14:11:33 +0000 Subject: [PATCH 4/9] Adding more timing statements --- time_integrators.f90 | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/time_integrators.f90 b/time_integrators.f90 index 1f966e9..f9ace5f 100755 --- a/time_integrators.f90 +++ b/time_integrators.f90 @@ -96,6 +96,7 @@ subroutine imex_rk(vtk_print, save_nusselt) 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) @@ -120,6 +121,8 @@ subroutine imex_rk(vtk_print, save_nusselt) Ti (:,it) = tmp_T uyi (:,it) = tmp_uy end do + finish = OMP_GET_WTIME() + write(*,*) " - stage 1 mid timing: ", finish-start, "(s)" ! Compute K2hat start = OMP_GET_WTIME() call calc_explicit(2) @@ -129,6 +132,7 @@ subroutine imex_rk(vtk_print, save_nusselt) !::::::::::: ! 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) @@ -153,6 +157,8 @@ subroutine imex_rk(vtk_print, save_nusselt) Ti (:,it) = tmp_T uyi (:,it) = tmp_uy end do + finish = OMP_GET_WTIME() + write(*,*) " - stage 2 mid timing: ", finish-start, "(s)" ! Compute K3hat start = OMP_GET_WTIME() call calc_explicit(3) @@ -162,6 +168,7 @@ subroutine imex_rk(vtk_print, save_nusselt) !::::::::::: ! 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) @@ -186,6 +193,8 @@ subroutine imex_rk(vtk_print, save_nusselt) Ti (:,it) = tmp_T uyi (:,it) = tmp_uy end do + finish = OMP_GET_WTIME() + write(*,*) " - stage 3 mid timing: ", finish-start, "(s)" ! Compute K4hat start = OMP_GET_WTIME() call calc_explicit(4) @@ -194,6 +203,7 @@ subroutine imex_rk(vtk_print, save_nusselt) ! 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,:)) + & @@ -217,6 +227,8 @@ subroutine imex_rk(vtk_print, save_nusselt) ux(:,it) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow! end if end do + finish = OMP_GET_WTIME() + write(*,*) " - update sols timing: ", finish-start, "(s)" if (time == t_final) then exit @@ -231,12 +243,15 @@ subroutine imex_rk(vtk_print, save_nusselt) end if end if + 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 + finish = OMP_GET_WTIME() + write(*,*) " - nusselt calc timing: ", finish-start, "(s)" finish_overall = OMP_GET_WTIME() write(*,*) "overall timing: ", finish_overall-start_overall, "(s)" From 6127cb8e3a627c00be4ceddd60c72b2059558574 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sun, 18 Apr 2021 23:03:51 +0000 Subject: [PATCH 5/9] NaN nusselts --- input.data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/input.data b/input.data index ca77d7c..d1caea6 100755 --- a/input.data +++ b/input.data @@ -2,6 +2,6 @@ 5.0e+03, 1, 1.1 5.0, 0.1 -1.0, 1.0 -1600, 1275, 1 +1600, 2500, 1 0, 0, 0 0, 0, 0 From f73bc6dbf7bd29f232d2ccdd57fd87902cf4792f Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Tue, 20 Apr 2021 13:52:39 +0000 Subject: [PATCH 6/9] 4x loop unrolling --- input.data | 2 +- time_integrators.f90 | 455 ++++++++++++++++++++++++------------------- 2 files changed, 257 insertions(+), 200 deletions(-) mode change 100755 => 100644 time_integrators.f90 diff --git a/input.data b/input.data index d1caea6..9711f31 100755 --- a/input.data +++ b/input.data @@ -2,6 +2,6 @@ 5.0e+03, 1, 1.1 5.0, 0.1 -1.0, 1.0 -1600, 2500, 1 +1600, 1276, 1 0, 0, 0 0, 0, 0 diff --git a/time_integrators.f90 b/time_integrators.f90 old mode 100755 new mode 100644 index f9ace5f..d0bb01d --- a/time_integrators.f90 +++ b/time_integrators.f90 @@ -40,15 +40,15 @@ subroutine imex_rk(vtk_print, save_nusselt) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 @@ -71,36 +71,36 @@ subroutine imex_rk(vtk_print, save_nusselt) ! Time integration do ! while (time < t_final) - start_overall = OMP_GET_WTIME() + 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 - 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 + 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) @@ -112,31 +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 - 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 + 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) @@ -148,28 +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 - 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 + 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 @@ -184,77 +184,77 @@ 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 - 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,:)) + & + 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 - finish = OMP_GET_WTIME() - write(*,*) " - update sols timing: ", finish-start, "(s)" + 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 + end if - !call update_dt - ! Don't write vtk for now. - wvtk = .false. - if (wvtk) then + !call update_dt + ! Don't write vtk for now. + wvtk = .false. + 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 - start = OMP_GET_WTIME() - ! 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 - 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 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 @@ -297,18 +297,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) @@ -326,7 +326,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) @@ -351,7 +351,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) @@ -360,27 +360,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) @@ -416,21 +416,21 @@ subroutine calc_explicit(stage) 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() @@ -438,10 +438,10 @@ subroutine calc_explicit(stage) 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)" @@ -450,26 +450,83 @@ subroutine calc_explicit(stage) 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 +do j = 1,Ny,4 + ! 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+1,:) + tnlphi = nlphi(j+1,:) + tT = Ti(j+1,:) + tux = uxi(j+1,:) + tuy = uyi(j+1,:) + tphi = phii(j+1,:) + 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+1,:) = tnlT + nlphi(j+1,:) = tnlphi + Ti(j+1,:) = tT + uxi(j+1,:) = tux + uyi(j+1,:) = tuy + phii(j+1,:) = tphi + ! Bring everything to physical space + tnlT = nlT(j+2,:) + tnlphi = nlphi(j+2,:) + tT = Ti(j+2,:) + tux = uxi(j+2,:) + tuy = uyi(j+2,:) + tphi = phii(j+2,:) + 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+2,:) = tnlT + nlphi(j+2,:) = tnlphi + Ti(j+2,:) = tT + uxi(j+2,:) = tux + uyi(j+2,:) = tuy + phii(j+2,:) = tphi + ! Bring everything to physical space + tnlT = nlT(j+3,:) + tnlphi = nlphi(j+3,:) + tT = Ti(j+3,:) + tux = uxi(j+3,:) + tuy = uyi(j+3,:) + tphi = phii(j+3,:) + 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+3,:) = tnlT + nlphi(j+3,:) = tnlphi + Ti(j+3,:) = tT + uxi(j+3,:) = tux + uyi(j+3,:) = tuy + phii(j+3,:) = tphi end do finish = OMP_GET_WTIME() write(*,*) " - - l3 timing: ", finish-start, "(s)" @@ -477,11 +534,11 @@ subroutine calc_explicit(stage) ! 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)" @@ -489,19 +546,19 @@ subroutine calc_explicit(stage) ! 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)" @@ -511,25 +568,25 @@ subroutine calc_explicit(stage) 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) @@ -566,11 +623,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 @@ -620,12 +677,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. @@ -644,38 +701,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)) From 4837549ca2443a6a4f3357c1f6954ef3547ebe70 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Tue, 20 Apr 2021 14:05:05 +0000 Subject: [PATCH 7/9] reset input --- input.data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/input.data b/input.data index 9711f31..ca77d7c 100755 --- a/input.data +++ b/input.data @@ -2,6 +2,6 @@ 5.0e+03, 1, 1.1 5.0, 0.1 -1.0, 1.0 -1600, 1276, 1 +1600, 1275, 1 0, 0, 0 0, 0, 0 From b15415e91597179042c33952ed8e5ddb199c006d Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Tue, 20 Apr 2021 14:34:27 +0000 Subject: [PATCH 8/9] reset after loop unrolling test --- time_integrators.f90 | 59 +------------------------------------------- 1 file changed, 1 insertion(+), 58 deletions(-) diff --git a/time_integrators.f90 b/time_integrators.f90 index d0bb01d..c2a7c2b 100644 --- a/time_integrators.f90 +++ b/time_integrators.f90 @@ -450,7 +450,7 @@ subroutine calc_explicit(stage) nlT = CI*nlT start = OMP_GET_WTIME() -do j = 1,Ny,4 +do j = 1,Ny ! Bring everything to physical space tnlT = nlT(j,:) tnlphi = nlphi(j,:) @@ -470,63 +470,6 @@ subroutine calc_explicit(stage) uxi(j,:) = tux uyi(j,:) = tuy phii(j,:) = tphi - ! Bring everything to physical space - tnlT = nlT(j+1,:) - tnlphi = nlphi(j+1,:) - tT = Ti(j+1,:) - tux = uxi(j+1,:) - tuy = uyi(j+1,:) - tphi = phii(j+1,:) - 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+1,:) = tnlT - nlphi(j+1,:) = tnlphi - Ti(j+1,:) = tT - uxi(j+1,:) = tux - uyi(j+1,:) = tuy - phii(j+1,:) = tphi - ! Bring everything to physical space - tnlT = nlT(j+2,:) - tnlphi = nlphi(j+2,:) - tT = Ti(j+2,:) - tux = uxi(j+2,:) - tuy = uyi(j+2,:) - tphi = phii(j+2,:) - 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+2,:) = tnlT - nlphi(j+2,:) = tnlphi - Ti(j+2,:) = tT - uxi(j+2,:) = tux - uyi(j+2,:) = tuy - phii(j+2,:) = tphi - ! Bring everything to physical space - tnlT = nlT(j+3,:) - tnlphi = nlphi(j+3,:) - tT = Ti(j+3,:) - tux = uxi(j+3,:) - tuy = uyi(j+3,:) - tphi = phii(j+3,:) - 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+3,:) = tnlT - nlphi(j+3,:) = tnlphi - Ti(j+3,:) = tT - uxi(j+3,:) = tux - uyi(j+3,:) = tuy - phii(j+3,:) = tphi end do finish = OMP_GET_WTIME() write(*,*) " - - l3 timing: ", finish-start, "(s)" From 32688c75bd2567d323bf833cc75fc413ba678203 Mon Sep 17 00:00:00 2001 From: Michael Neuder Date: Fri, 23 Apr 2021 11:10:25 -0400 Subject: [PATCH 9/9] Vtk experiments --- input.data | 4 ++-- time_integrators.f90 | 2 -- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/input.data b/input.data index ca77d7c..894b925 100755 --- a/input.data +++ b/input.data @@ -2,6 +2,6 @@ 5.0e+03, 1, 1.1 5.0, 0.1 -1.0, 1.0 -1600, 1275, 1 -0, 0, 0 +64, 51, 1 0, 0, 0 +1, 0, 0 diff --git a/time_integrators.f90 b/time_integrators.f90 index c2a7c2b..8149103 100644 --- a/time_integrators.f90 +++ b/time_integrators.f90 @@ -235,8 +235,6 @@ subroutine imex_rk(vtk_print, save_nusselt) end if !call update_dt - ! Don't write vtk for now. - wvtk = .false. if (wvtk) then if (mod(nti,vtk_print) == 0) then call write_to_vtk(nti, .false.) ! false = Fourier space