diff --git a/src/common.f90 b/src/common.f90 index 31f5f29c..9d0c5a86 100644 --- a/src/common.f90 +++ b/src/common.f90 @@ -15,6 +15,7 @@ module m_common integer :: n_groups_x, n_groups_y, n_groups_z real(dp) :: Lx, Ly, Lz real(dp) :: dx, dy, dz + real(dp) :: nu, dt integer :: nproc_x = 1, nproc_y = 1, nproc_z = 1 character(len=20) :: BC_x_s, BC_x_e, BC_y_s, BC_y_e, BC_z_s, BC_z_e integer :: poisson_solver_type diff --git a/src/solver.f90 b/src/solver.f90 index b09fab24..825c1481 100644 --- a/src/solver.f90 +++ b/src/solver.f90 @@ -84,8 +84,8 @@ function init(backend, time_integrator, xdirps, ydirps, zdirps, globs) & real(dp), allocatable, dimension(:, :, :) :: u_init, v_init, w_init integer :: dims(3) - real(dp) :: dx, dy, dz - integer :: nx, ny, nz + real(dp) :: x, y, z + integer :: nx, ny, nz, i, j, b, ka, kb, ix, iy, iz, sz solver%backend => backend solver%time_integrator => time_integrator @@ -108,6 +108,30 @@ function init(backend, time_integrator, xdirps, ydirps, zdirps, globs) & v_init = 0 w_init = 0 + solver%dt = globs%dt + solver%backend%nu = globs%nu + + sz = dims(1) + nx = globs%nx_loc; ny = globs%ny_loc; nz = globs%nz_loc + do ka = 1, nz + do kb = 1, ny/sz + do j = 1, nx + do i = 1, sz + ! Mapping to ix, iy, iz depends on global group numbering + ix = j; iy = (kb - 1)*sz + i; iz = ka + x = (ix - 1)*globs%dx + y = (iy - 1)*globs%dy + z = (iz - 1)*globs%dz + + b = ka + (kb - 1)*xdirps%n + u_init(i, j, b) = sin(x)*cos(y)*cos(z) + v_init(i, j, b) =-cos(x)*sin(y)*cos(z) + w_init(i, j, b) = 0 + end do + end do + end do + end do + call solver%backend%set_fields( & solver%u, solver%v, solver%w, u_init, v_init, w_init & ) @@ -115,13 +139,10 @@ function init(backend, time_integrator, xdirps, ydirps, zdirps, globs) & deallocate(u_init, v_init, w_init) print*, 'initial conditions are set' - nx = globs%nx_loc; ny = globs%ny_loc; nz = globs%nz_loc - dx = globs%dx; dy = globs%dy; dz = globs%dz - ! Allocate and set the tdsops - call allocate_tdsops(solver%xdirps, nx, dx, solver%backend) - call allocate_tdsops(solver%ydirps, ny, dy, solver%backend) - call allocate_tdsops(solver%zdirps, nz, dz, solver%backend) + call allocate_tdsops(solver%xdirps, nx, globs%dx, solver%backend) + call allocate_tdsops(solver%ydirps, ny, globs%dy, solver%backend) + call allocate_tdsops(solver%zdirps, nz, globs%dz, solver%backend) select case (globs%poisson_solver_type) case (POISSON_SOLVER_FFT) diff --git a/src/xcompact.f90 b/src/xcompact.f90 index 1600d4ad..3a2726fa 100644 --- a/src/xcompact.f90 +++ b/src/xcompact.f90 @@ -61,6 +61,9 @@ program xcompact ! read ns from the input file globs%nx = 512; globs%ny = 512; globs%nz = 512 + globs%dt = 0.001_dp + globs%nu = 1._dp/1600._dp + ! set nprocs based on run time arguments globs%nproc_x = 1; globs%nproc_y = 1; globs%nproc_z = 1 @@ -123,8 +126,6 @@ program xcompact print*, 'OpenMP backend instantiated' #endif - backend%nu = 1._dp - allocate(u(SZ, globs%nx_loc, globs%n_groups_x)) allocate(v(SZ, globs%nx_loc, globs%n_groups_x)) allocate(w(SZ, globs%nx_loc, globs%n_groups_x))