diff --git a/vlaplex.f90 b/vlaplex.f90 index d82a49d..3178b5f 100644 --- a/vlaplex.f90 +++ b/vlaplex.f90 @@ -59,7 +59,7 @@ program VlaPlEx character(:), allocatable:: bc_file real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old - real(dp), allocatable, dimension(:,:):: f0 ! Boundary at r = x_0 + real(dp), allocatable, dimension(:):: f0 ! Boundary at r = x_0 real(dp), allocatable, dimension(:,:):: n_i real(dp), allocatable, dimension(:,:):: u_i real(dp), allocatable, dimension(:):: E_i @@ -158,13 +158,13 @@ program VlaPlEx nz = 2 ! Allocate vectors - allocate(f_i(1:nr,1:nv,1:nz), f_i_old(1:nr,1:nv,1:nz)) - allocate(n_i(1:nr,1:nz)) - allocate(u_i(1:nr,1:nz), E_i(1:nr), T_i(1:nr,1:nz)) + allocate(f_i(1:nz,1:nr,1:nv), f_i_old(1:nz,1:nr,1:nv)) + allocate(n_i(1:nz,1:nr)) + allocate(u_i(1:nz,1:nr), E_i(1:nr), T_i(1:nz,1:nr)) allocate(Zave(1:nr)) allocate(n_e(1:nr)) allocate(phi(1:nr), phi_old(1:nr), E(1:nr)) - allocate(fCum_i(1:nv,1:nz)) + allocate(fCum_i(1:nz,1:nv)) f_i = 0.0_dp f_i_old = 0.0_dp n_i = 0.0_dp @@ -217,7 +217,7 @@ program VlaPlEx phi0 = 1.0e2_dp / phi_ref ! Dirichlet phi(1) = phi0 ! Dirichlet ! phi0 = phi(1) ! Neumann - allocate(f0(j0:nv,1:nz)) + allocate(f0(j0:nv)) f0 = 0.0_dp ! Output initial values @@ -228,7 +228,7 @@ program VlaPlEx ! call writeOutputF(t, dt, nr, r, nv, v, f_i_old) call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i) call writeOutputPhi(t, dt, nr, r, phi, E, n_e) - call writeOutputMom(t, dt, nr, r, n_i(:,1), u_i(:,1), T_i(:,1), Zave) + call writeOutputMom(t, dt, nr, r, n_i(1,:), u_i(1,:), T_i(1,:), Zave) ! Main loop do t = 1, nt @@ -238,14 +238,14 @@ program VlaPlEx u_bc = sqrt(Zave_bc * Temp_bc) do iz = 1, nz ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) - f0(j0:nv,iz) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) - f0(:,iz) = f0(:,iz) * n_bc / (sum(f0(:,iz))*dv) + f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) + f0 = f0 * n_bc / (sum(f0)*dv) ! Boundary conditions ! r = r0, v>0 - f_i_old(1,j0:nv,iz) = f0(:,iz) - f_i(1,j0:nv,iz) = f_i_old(1,j0:nv,iz) - T_i(1,iz) = Temp_bc + f_i_old(iz,1,j0:nv) = f0 + f_i(iz,1,j0:nv) = f_i_old(iz,1,j0:nv) + T_i(iz,1) = Temp_bc end do T_e = Temp_bc print *, 'Time: ', time * t_ref @@ -258,36 +258,36 @@ program VlaPlEx Zave(1) = Zave_bc ! r = rf, v<0 - f_i_old(nr,1:j0-1,:) = 0.0_dp - f_i(nr,1:j0-1,:) = f_i_old(nr,1:j0-1,:) + f_i_old(:,nr,1:j0-1) = 0.0_dp + f_i(:,nr,1:j0-1) = f_i_old(:,nr,1:j0-1) ! set edge velocities to 0 - f_i_old(:,1,:) = 0.0_dp - f_i_old(:,nv,:) = 0.0_dp + f_i_old(:,:,1) = 0.0_dp + f_i_old(:,:,nv) = 0.0_dp ! Advect in the r direction !$omp parallel do do iz = 1, nz do i = 1, nr ! Advect negative velocity if (i < nr) then - f_i(i,1:j0-1,iz) = f_i_old(i,1:j0-1,iz) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(i+1,1:j0-1,iz) - & - r(i )**2*f_i_old(i ,1:j0-1,iz)) + f_i(iz,i,1:j0-1) = f_i_old(iz,i,1:j0-1) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(iz,i+1,1:j0-1) - & + r(i )**2*f_i_old(iz,i ,1:j0-1)) end if ! Advect positive velocity if (i > 1) then - f_i(i,j0:nv,iz) = f_i_old(i, j0:nv,iz) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(i , j0:nv,iz) - & - r(i-1)**2*f_i_old(i-1, j0:nv,iz)) + f_i(iz,i,j0:nv) = f_i_old(iz,i, j0:nv) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(iz,i , j0:nv) - & + r(i-1)**2*f_i_old(iz,i-1, j0:nv)) end if - n_i(i,iz) = sum(f_i(i,:,iz))*dv - if (n_i(i,1) > 1.0e-10_dp) then - u_i(i,iz) = sum(v(:) *f_i(i,:,iz))*dv / n_i(i,iz) - E_i(i) = sum(v(:)**2*f_i(i,:,iz))*dv / n_i(i,iz) - T_i(i,iz) = 2.0_dp*E_i(i) - 2.0_dp*u_i(i,iz)**2 + n_i(iz,i) = sum(f_i(iz,i,:))*dv + if (n_i(1,i) > 1.0e-10_dp) then + u_i(iz,i) = sum(v(:) *f_i(iz,i,:))*dv / n_i(iz,i) + E_i(i) = sum(v(:)**2*f_i(iz,i,:))*dv / n_i(iz,i) + T_i(iz,i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(iz,i)**2 Zave(i) = Zave_bc else - u_i(i,iz) = 0.0_dp - T_i(i,iz) = 0.0_dp + u_i(iz,i) = 0.0_dp + T_i(iz,i) = 0.0_dp Zave(i) = 0.0_dp end if @@ -296,7 +296,7 @@ program VlaPlEx !$omp end parallel do ! Assume quasi-neutrality to start iterating - n_e = Zave * n_i(:,1) + n_e = Zave * n_i(1,:) !do iz = 1, nz ! n_e = n_e + Zave * n_i(:,iz) !end do @@ -322,7 +322,7 @@ program VlaPlEx !do iz = 1, nz ! b = b - (Zave * n_i(:,iz)) !end do - b = -(Zave * n_i(:,1) - n_e) + b = -(Zave * n_i(1,:) - n_e) ! Apply boundary conditions b(1) = phi0 ! Dirichlet ! b(nr) = 0.0_dp ! Dirichlet @@ -383,18 +383,18 @@ program VlaPlEx ! i = 1, v<0 i = 1 if (E(i) >= 0.0_dp) then - f_i(i,2:j0-2,iz) = f_i_old(i,2:j0-2,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:j0-2,iz) - f_i_old(i,1:j0-3,iz)) + f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,2:j0-2) - f_i_old(iz,i,1:j0-3)) else - f_i(i,2:j0-2,iz) = f_i_old(i,2:j0-2,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:j0-1,iz) - f_i_old(i,2:j0-2,iz)) + f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,3:j0-1) - f_i_old(iz,i,2:j0-2)) end if ! i = 2, nr-1; all v !$omp parallel do do i = 2, nr-1 if (E(i) >= 0.0_dp) then - f_i(i,2:nv-1,iz) = f_i_old(i,2:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:nv-1,iz) - f_i_old(i,1:nv-2,iz)) + f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,2:nv-1) - f_i_old(iz,i,1:nv-2)) else - f_i(i,2:nv-1,iz) = f_i_old(i,2:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:nv,iz) - f_i_old(i,2:nv-1,iz)) + f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,3:nv) - f_i_old(iz,i,2:nv-1)) end if @@ -403,9 +403,9 @@ program VlaPlEx ! i = nr, v>=0 i = nr if (E(i) >= 0.0_dp) then - f_i(i,j0+1:nv-1,iz) = f_i_old(i,j0+1:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+1:nv-1,iz) - f_i_old(i,j0:nv-2,iz)) + f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,j0+1:nv-1) - f_i_old(iz,i,j0:nv-2)) else - f_i(i,j0+1:nv-1,iz) = f_i_old(i,j0+1:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+2:nv,iz) - f_i_old(i,j0+1:nv-1,iz)) + f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,j0+2:nv) - f_i_old(iz,i,j0+1:nv-1)) end if end do @@ -413,14 +413,14 @@ program VlaPlEx ! Reset values for next iteration f_i_old = f_i do iz = 1, nz - fCum_i(:,iz) = fCum_i(:,iz) + f_i_old(rCum_index,:,iz) + fCum_i(iz,:) = fCum_i(iz,:) + f_i_old(iz,rCum_index,:) end do ! Write output if (mod(t,everyOutput) == 0 .or. t == nt) then ! call writeOutputF(t, dt, nr, r, nv, v, f_i_old) call writeOutputPhi(t, dt, nr, r, phi, E, n_e) - call writeOutputMom(t, dt, nr, r, n_i(:,1), u_i(:,1), T_i(:,1), Zave) - call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i(:,1)) + call writeOutputMom(t, dt, nr, r, n_i(1,:), u_i(1,:), T_i(1,:), Zave) + call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i(1,:)) end if